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We explore application of the domain wall fermion formalism of lattice QCD to calculate the 
K — > TTiv decay amplitudes in terms of the K + — > tt + and K° — > hadronic matrix elements 
through relations derived in chiral perturbation theory. Numerical simulations are carried out in 
quenched QCD using domain- wall fermion action for quarks and an RG-improved gauge action for 
giuons on a 16 3 x 32 x 16 and 24 3 x 32 x 16 lattice at j3 = 2.6 corresponding to the lattice spacing 
1/a ~ 2 GeV. Quark loop contractions which appear in Penguin diagrams are calculated by the 
random noise method, and the A7 =1/2 matrix elements which require subtractions with the quark 
loop contractions are obtained with a statistical accuracy of about 10%. We investigate the chiral 
properties required of the K + — > ix + matrix elements. Matching the lattice matrix elements to 
those in the continuum at fi = 1/a using the perturbative renormalization factor to one loop order, 
and running to the scale /j, = m c = 1.3 GeV with the renormalization group for Nf = 3 flavors, we 
calculate all the matrix elements needed for the decay amplitudes. 

With these matrix elements, the AI = 3/2 decay amplitude Re^2 shows a good agreement with 
experiment after an extrapolation to the chiral limit. The A7 = 1/2 amplitude Re^lo, on the other 
hand, is about 50-60% of the experimental one even after chiral extrapolation. In view of the 
insufficient enhancement of the AI = 1/2 contribution, we employ the experimental values for the 
real parts of the decay amplitudes in our calculation of e' /e. The central values of our result indicate 
that the A7 = 3/2 contribution is larger than the AI = 1/2 contribution so that e'/e is negative 
and has a magnitude of order 10 -4 . We discuss in detail possible systematic uncertainties, which 
seem too large for a definite conclusion on the value of e'/e. 

PACS numbers: 11.15Ha, 12.38Gc 



I. INTRODUCTION 

Understanding non-leptonic weak processes of kaon, in particular the K — > 7T7t decay, represents one of the keys to 
establishing the Standard Model and probing the physics beyond it. This decay exhibits two significant phenomena, 
namely, the AI = 1/2 rule, which is a large enhancement of the decay mode with AI =1/2 relative to that with 
AI = 3/2, and the direct CP violation 0,0 > which is naturally built in the Model for three or more families of 
quarks y|. While both of these phenomena are well established by experiments, theoretical calculations with sufficient 
reliability that allow examinations of the Standard Model predictions against the experimental results are yet to be 
made. The main reason for this status is the difficulty in calculating the hadronic matrix elements of local operators 
which appear in the effective weak Hamiltonian for the decay amplitudes. At the energy scales relevant for these 
operators, analytic treatments such as the l/N c expansion are not sufficiently powerful to reliably evaluate the effect 
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of the strong interactions in the matrix elements. In fact, the AI =1/2 rule, which is supposed to arise from QCD 
effects, has not been quantitatively explained by analytic methods so far. With these backgrounds, Monte Carlo 
simulations of lattice QCD provide a hopeful method for the calculation of the decay amplitudes. 

A natural framework for theoretical calculations of the decay amplitudes is provided by the effective weak Hamil- 
tonian, Hw, which follows from an operator product expansion (OPE) of weak currents Q: 

H W = ^KsKTd E Wi(jM)Qi(jj). (LI) 

Here the Wilson coefficients Wi contain effects of the energy scales higher than /i so that they can be calculated 
perturbatively. Non-perturbative QCD effects are contained in the matrix elements of the local operators Qi, and 
the calculation of these matrix elements, often called hadronic matrix elements (HME), is the task of lattice QCD 
0- IE 0- • Our a i m m this paper is to report on our attempt to obtain these matrix elements through numerical 
simulations of lattice QCD using the domain wall formalism I! EH EH for quarks. 

The amplitudes for K — > tttt decay with AI =1/2 and 3/2 are written as the matrix elements of Hyy, 

((^) I \H w \K°) = A I e lS >, (1.2) 

where the subscript 7 = or 2 denotes the isospin of the final state corresponding to AJ = 1/2 or 3/2, and 5i is the 
phase shift from final state interactions tttt — > tttt caused by QCD effects. The AI = 1/2 rule, which is one of the 
focuses of our calculation, is described by the ratio of isospin amplitudes Ai\ 

u -i = ^£ « 22.2. (1.3) 
ReA 2 

Another focus is the parameter e'/e of direct CP violation in the Standard Model. The recent experimental results 
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_ J (20.7 ± 2.8) • 10~ 4 (KTeV) [I] , . 

(15.3 ±2.6) • 10" 4 (NA48) [2] ' v > 



In the numerical simulation of lattice QCD, matrix elements are generally extracted from Euclidean correlation 
functions of the relevant operators and those which create the initial and final states in their lowest energy levels. 
For sufficiently large Euclidean time distances, excited states damp out and the matrix elements of the lowest energy 
states are left. In fact, the kaon B-parameter Bk has been successfully obtained from the three-point correlation 
function of K° and K and an insertion of the AS = 2 weak Hamiltonian E3- However, in the calculation of 
the four-point function, (7r(^)7r(ti ) Hw(tn )K(tx)), necessary for the K — > tttt decay, there is a severe limitation as 
pointed out by Maiani and Testa They have shown that it is difficult to obtain the matrix elements unless the 
momentum of each of the two pions in the final state is set to zero. 

One of the ways to overcome the difficulty pursued in the past is to calculate the matrix elements with the two pions 
at rest, allowing a nonzero energy transfer AE = 2771^ — ttik at the weak operator. This generally causes mixings of 
unphysical lower dimension operators through renormalization, which has to be removed. (See Ref. Q and references 
therein.) Furthermore, the unphysical amplitudes obtained with AE ^ need to be extrapolated to physical ones 
by use of some effective theories such as chiral perturbation theory. Due to these problems and numerical difficulties 
of extracting reasonable signals from four-point functions, this approach has not been successful for the A J = 1/2 
amplitude despite many efforts over the years E3- F° r the AI = 3/2 amplitude for which the operator mixing is 
absent, on the other hand, a recent study has obtained a result in agreement with experiment [Tr| . 

Several prop osals have been presented over the years for extracting the physical amplitude from the four-point 
functions |l7tll8l.ll9|. Feasibility studies for implementing them in practical simulations are yet to come, however. 

In this paper we explore a method proposed by Bernard et al. |l5j | which is alternative to calculating the three-point 
function. In this method, which we shall call as reduction method, chiral perturbation theory (%PT) is used to relate 
the matrix elements for K — > tttt to those for K — > tt and K — > (vacuum) , and the latter amplitudes are calculated in 
lattice QCD. Since this calculation involves only three- and two-point correlation functions, the Maiani- Testa problem 
mentioned above is avoided. Statistical fluctuations are also expected to be diminished compared with the case of 
four-point correlation functions. 

Early attempts with this method |l'l| encountered large statistical fluctuations in the correlation functions so that 
meaningful results were difficult to obtain. For the Wilson fermion action or its 0{a) improved version, there is an 
added difficulty that the mixing of operators of wrong chirality caused by explicit chiral symmetry breakin g o f the 
action has to be removed. The mixing problem has been resolved only for the AI = 3/2 operators so far |20l I2H |22 | ■ 
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The first results on the AI =1/2 rule and e' /e calculated with this method were recently reported |2^| using the 
staggered fermion action which keeps t/(l) subgroup of chiral symmetry. In this work, however, a large dependence 
of the AI = 3/2 amplitude on the meson mass was seen, which made the chiral extrapolation difficult. Moreover, 
large uncertainties due to perturbative renormalization factors depending on the value of the matching point were 
reported. Hence clear statements on the viability of the method were difficult to make from this work. 

In this article we report on our attempt to apply the domain- wall fermion formalism of lattice QCD to the 

calculation of K — > mr decay amplitudes in the context of the reduction method. A major advantage of this approach 
over the conventional fermion formalisms is that full chiral symmetry can be expected to be realized for sufficiently 
large lattice sizes in the fifth dimension. Good chiral property of one of the K — > 7r matrix elements, equivalent to the 
kaon B parameter, was observed in the pioneering application of the formalism j24[. Detailed investigations into the 
realization of the chiral limit have been made in the quenched approximation for the plaquette and a renormalization 
group (RG) improved gluon actions l26l |27|. It was found that the use of RG- improved action leads to much 
better chiral properties compared to the case of the plaquette action for similar lattice spacings [2(j. This prompts 
us to adopt the RG-improved action in our simulation. 

Another possible advantage of the domain wall formalism is 0(a 2 ) scaling violation from the fermion sector as 
opposed to 0(a) for the Wilson case. Indeed our domain wall fermion calculation of Bk |28( exhibits only small 
scaling violation. The magnitude of violation is much smaller compared to the staggered fermion case [2^ which is 
also expected to be 0(a 2 ). An improved scaling behavior may be enhanced with the use of the RG-improved gluon 
action. 

This article is organized as follows. In Section [HJ we summarize the main points of the %PT reduction method. 
For the construction of the formulae which relate the matrix elements for K — » 7r and the K — > tttt decay amplitudes, 
the relations between the four quark operators Qi and xPT operators are considered at tree level on the basis of 
chiral transformation properties. The necessity of chiral symmetry on the lattice is emphasized. In Section IIIII we 
summarize the details of our numerical simulation procedure. We discuss the form of lattice actions and the choice 
of optimal set of simulation parameters from the point of view of chiral properties. Some of the technical issues are 
also explained including renormalization of the four-quark operators and RG-running of the matrix elements to the 
relevant energy scale. The numerical results are reported in Section llVl and Ivl The former contains results of hadronic 
matrix elements. In particular, we show that the subset of K — > tt matrix elements which are expected to vanish in 
the chiral limit satisfy this requirement. We then present the physical matrix elements and combine them with the 
Wilson coefficients, which are already calculated perturbatively. This leads us to results for the AI =1/2 rule and 
e'/e. Our conclusions are given in Section IVT1 

A preliminary report of the present work was presented in Ref. [3(J. We refer to Refs. [3lll32]| for a similar attempt, 
and Refs. 8, 33 for reviews. 



II. CHIRAL PERTURBATION THEORY REDUCTION METHOD 



A. LOCAL OPERATORS 

We carry out our analyses choosing the energy scale /i in the OPE for the weak Hamiltonian (|I.1|) equal to the charm 
quark mass m c = 1.3 GeV. In this case only u, d and s quarks appear in the local four-quark operators. Conventionally 
these operators are written as 



Qi = 




-f5)Ub][u b J^(l - -f5)d a ], 


(III) 


Q2 = 




IhhaWubl^l ~ lb)db], 


(II.2) 


Qz = 


[s 7 M (l - 


l5)d a ]^2,[qbl^ ~ 75)%], 


(II.3) 


Qi = 




75)4] ^2\qbi^{i - ih)q a ], 
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(II.4) 


Q 5 = 


[s a 7 M (l - 
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(II.5) 


Q 6 = 


[s Q 7 M (l - 


75)4] ^[<ftOv(l + 75)?aL 
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(II.6) 


Q7 = 




- 7sR] ^2 e 9[9b7^( 1 + 75)96], 


(II.7) 



9 



4 



Qs = ^[sal^l-l^d^^eqlq^^l+^qa], (II.8) 

8 

Qg = ^[v^l ~ "f5)da\ ^ejg^l ~ 75)%], (II.9) 

9 

QlO = ■^[saj f ,(l~j5)db}^2e q [q b j f _ l (l--f 5 )q a }, (11.10) 

« 

where the indices a, 6 denote color, and the summation over q appearing in Q 3 to Q10 runs over the three light flavors, 
q = u,d,s, with the charge e u = 2/3 and e<j = e s — —1/3. 

With the use of Fierz rearrangements, one can derive the relations, 

Qi = Q2 + Q3-Q1, (Il.ll) 

Qo = (11.12) 

Qio = \Qi-\q± = Qi - \q 5 + \Qi- (11.13) 

Hence Q4, Qg, Q10 are not independent operators. We emphasize that these relations do not hold in general d- 
dimensions where Fierz rearrangements cannot be used. 

In terms of the irreducible representations of chiral SU(3)l <8> SU(3) r group, Qi's are classified as 

(27l, 1 B ) © (8i, 1*), (11.14) 

$l,1r), (H.15) 

(8i,8 fl ). (11.16) 



Qi, Q2, Q^i Qw 
Q3, Q4, Q5, Q& 



(8lAr) 




= (sd) L (uu) L - {su) L (ud) L , 




(8lAr) 


x 2 


= {sd) L [{uu) L + 2{dd) L + 2(ss) L ] + {su) L 


{ud) L 


(27i,l fl ) 


x 3 


= (sd) L [2{uu) L - (dd) L - (ss)l\ + 2{su) L 


{ud)h 


(8lAr) 


Yx 


= (sd) L [(uu) R + (dd) R + {ss) R ] , Yi c , 




(8l,8 k ) 


Y 2 


= (sd) L [2{uu) R - (dd) R - {ss) R ] , Y 2 C , 





The operators Qi(i = 1, • • • , 10) are invariant under CPS symmetry, i.e., the product of CP transformation and d <-» s 
interchange. A basis of operators which are irreducible under chiral symmetry and invariant under CPS is given by 

(11.17) 

(11.18) 

(11.19) 

(11.20) 

(11.21) 

where (sd)^ = sj^(l — 7s)d and {sd) R = sj^(l +j$)d. The color and spinor indices are summed within each current 
except for Yf for which the color summation is taken across the two currents. While X^s have the Lorentz structure 
of L <£> L, Yi's have that of L ® R. All the independent local operators are written as linear combinations of these 
operators : 

Qi = l -X 1 + ^X 2 + l -X^ (11.22) 
Q2 = - l -X 1 + ^X 2 + l -X z , (11.23) 

Qs = \x 1+ l -X 2: (11.24) 

Q 5 = Y u (11.25) 
Q 6 = Y x \ (11.26) 

Qr = \y 2 , (11.27) 
Qs = \y 2 c - (11.28) 



The expressions for the dependent operators (3.4.9,10 are easily derived using i|II.ll|) -( |II.13|) . 

The final states in the K — > irw decay can have either isospin / = or 2, i.e. AI = 1/2 or 3/2. Hence Qi's are 
decomposed as 

Q 1 = qT ) + Q ( P- (H.29) 

This decomposition is accomplished by constructing another basis of irreducible representations with the intrinsic 
isospin /. The details are described in Appendix lAl 
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B. CHIRAL PERTURBATION THEORY 

In the low energy region of strong interactions, the octet of pseudo scalar mesons 7r°, ir^, K°, ~K° , K ± , rj play a 
principal role as the Nambu-Goldstone bosons of spontaneously broken chiral symmetry SU(3)l®SU(3)h — ► SU(3)y. 
In chiral perturbation theory (xPT) as a low energy effective theory of QCD, these Nambu-Goldstone boson fields 
are used to parametrize the broken axial symmetry, and we collect them in a 3 x 3 matrix, 

E = (e**/') , 

a 

where A a are Gell-Mann matrices, and / is the decay constant. Under SU(3)l <8> SU(3)r chiral transformation, 
E € SU (3) transforms as 

E -» g R Vgl , Et -» g L Y^ g R . (11.32) 
The chiral Lagrangean to the lowest order, with the additional mass term, is given by 

C x = ^tr^E^E) - ^tr[M(Et + E)], (11.33) 

where M = (2 Bo) ■ diag[m„, m<j, m s ] denotes the quark mass matrix and Bq is a parameter. In terms of E, the left- 
and right-handed currents are given by 

M = \f{d^ ■ E)| , (R,)) = l -f(d^ ■ Et)|, (11.34) 

respectively. 

The idea of xPT reduction method by Bernard et al. [w| is to relate the hadronic matrix elements for K — > tttt 
decays to those for K — > tt and K ^ (vacuum) using xPT, and calculate the latter through numerical simulations of 
lattice QCD. As the first step of xPT reduction method, we construct operators in %PT which correspond to X^s and 
YiS in QCD i.e., those which transform under the same irreducible representations of SU(3)l® SU(3)r and invariant 
under CPS symmetry. In the following, we discuss the case of {(27l, 1r), (8l, and (8^,8^) representations 
separately. 



K~ 



7T+ K+ 



(11.30) 
(11.31) 



C. REDUCTION METHOD FOR (27 l ,1r) AND (8 l ,1r) OPERATORS 

For the irreducible representations (27 l, 1r) and (8l, 1r), which cover Q\, . . . , Qq, Qg and Qio, the product of left 
handed currents (Z/ M )*(L M )f is one of the candidates for the operator to the lowest order in xPT. An explicit form of 
the operators, which are also CPS invariant, is given by 

(8 l ,1r) ■ A (11.35) 
(27 L ,1 R ) : C = 3(L M )§(£ p )l + 2(L |i )J(£ M )?. (11.36) 

where A corresponds to X\ or X2, while C is the counterpart of X3. The latter is decomposed into two parts with 
1 = and 2 in the same way as X 3 (see Appendix lA"j l : 

C = Ic(°' + |c( 2 ), (11.37) 

where 

C (0) = (L^KL^l + (^1(^)1 + 2(^)1(^)1-3(^)1(^)1 (11.38) 

= (L,)\(L,)l + (L,)\(L,)l-(L,)l(L,)l (11.39) 

In addition to the operators above, there is another (8l, 1.r) operator which is allowed from CPS invariance: 

(8i,l R ) : B = (EM + MEt) 2 , 
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B a (m s + m d ) (S + E+)g - B (m s - m d ) (S - E f 



m d p 6 



(11.40) 



where = + Ly)j1 and A M = (i?^ — L^)/2 are vector and axial vector currents with and i? M defined in 



(|II.34|I . The equation of motion for £ is used to derive the third line from the second line in l|II.40(l 

The counterpart of this operator for QCD can be obtained easily by SU(3)l ® SU(3)r and CPS symmetry, 



;sub 



= d, 



(m s + m d )sd - (m s - m d )sj 5 d 
m., + m d _ , m.„ - m d 



md _ m 

S7 M d 

m d m 



m d 



(11.41) 



where the equation of motion for s and d quark fields is used. 

For physical K — > 7T7t processes, Q S ub, and hence B, does not contribute since these operators are total derivative 
of local operators and the energy-momentum injected at the weak operator vanishes. However, for the unphysical 
processes such as K ~ > 7r and K — > (vacuum) which we are to calculate on the lattice, the matrix elements of Q su b 
or B does not vanish due to a finite energy-momentum transfer for m s ^ m d . Therefore, a mixing between Q^s and 
Q su b in K — > 7r matrix elements exists which should be removed. We should also note that this mixing inevitably 
arises in the case of m d = m s , as is often chosen in numerical simulations on the lattice, since Q S ub is not a total 
divergence for this case. 

We assume that there are linear relations in the sense of matrix elements between the local operators {Qi (i — 
1, . . . , 6, 9, 10), Q S ub} and {A, B, C} which belong to the same representations, i.e. {(27l, 1r), (8l, 1r)} : 



Q S ub = rB, 

q| 2) = c | 2) c( 2 ), 



c<°> 



(11.42) 
(11.43) 
(11.44) 



0, if H 



where the coefficients ai,bi,c\ and r are unknown parameters. Taking the matrix elements of the two sides of (|II.42ll . 
and (ITTT441 for K° 



Q 



(0) 



Q 



(0) 



■ 7T+ and K° 
- UiQ sub K° 

>(2) 



one obtains 



0. 

2pk ■ P-n 



Q 



(0) 



(2) 



P 

2PK ■ Pit (2 



c ! 0) ) + o(/), 
+ o(p 4 ), 

.7 " 

3 (m^-m2)(a l - C f ,) ) + 0(p 4 ), 



72 
/ 

72 
/ 



G(P 4 



(11.45) 
(11.46) 
(11.47) 

(11.48) 

(11.49) 



where = fe^/r in (|II.45|) and (|II.46() . pk and are the momenta of kaon and pion, respectively, and p denotes 



either of them. In l|II.48|) and l|II.49|) . mx and are the physical meson masses. After eliminating — cf ^ from 



(|II.46|I and (|II.48() . we arrive at the relation between 



Q 



(0) 



K° ) and ( 7T 



Q 



(0) 



K + ) in the 1 = case: 



(0) 



(m 



V2f( P 



K - Pit) 



, sub 



(Xj 





Ql 0) 






Qsub 


K0) 



6, 9,10. 



(11.50) 
(11.51) 



The K — > (vacuum) matrix elements are used only to determine the oli 's which govern the subtraction of unphysical 
contributions originating from Q S ub- The relation for the 1 = 2 case is derived in the same way from l|II.47j) and 
(HL49ll : 



^-\qP\k°) = 



( m K ~ m l) 
V2f(PK ■ Ptv) 



n+\Q^\K+)+0(p 2 ), 



i = 6,9, 10. 



(11.52) 
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Let us note that the essential point of the reduction method is a calculation of the parameters a, — c| ' and cf^ 
from K — > 7r three-point correlation functions in numerical simulations of lattice QCD. Since these parameters appear 
in (|II.46|) and (|II.47|) as the coefficients of pk • Pit, their values are sensitive to the chiral properties of the K — > it 
matrix elements on the left hand side of these equations. Hence SU(3)l <& SU(3)r chiral symmetry on the lattice is 
an indispensable requirement for a successful calculation using this method. 

D. REDUCTION METHOD FOR (8 L ,8 R ) OPERATORS 

In order to construct (&l,8r) operators in xPT, we observe that (E)*-(£t)^ transforms as (8r, 8l) where 
(j, k) and correspond to 8l and 8r, respectively. One finds a CPS invariant operator 



D = 3£i(£t)f 

as the counterpart of Y2 ■ The decomposition into the 1 = and 2 part is given by 

D = D(°»+D< 2 ), 
D<°) = 2£i(£t)2_ s 2 (s t ) i + s 2 (Et) 2 j 

D (2) = £j(£t)? + £l(£t)l-£l(£t)2. 
Assuming linear relations between {Q^\ Q^} and D^'s, 

Ql 0) =dfr>M , Q^ = 4^HW (i = 7,8), 
with the unknown parameters d\ s, we take the matrix elements of the two sides for K — » irir and K — > ir to obtain 



(tt+ |Qf j K+) = 4di° V/ 2 + Ob 2 ) , #°) = -2V2d 4 (0) // 3 + G(P 2 ) 

q| 2) k°U-V24 2) // 3 + o(p 2 ) 



7(2) If2 



These relations lead to the reduction formulae for (8^,8^) operators, namely, 



Q 



a) 



K~ 



0(p 2 ), 



i = 7, 



which is common for the 1 = and 2 components. 



III. DETAILS OF SIMULATIONS 
A. LATTICE ACTIONS 



(11.53) 

(11.54) 
(11.55) 
(11.56) 

(11.57) 



(11.58) 
(11.59) 

(11.60) 



The RG-improved gauge action we use is defined by 

.5 (i 



1 plaqucttc 



pi + ci Yl 



rtg 



1x2 rectangle 



(III.1) 



where the coefficients of the plaquette and 1x2 Wilson loop terms take the values Cq = 3.648 and C\ = —0.331 (37J. 
This action is expected to lead to a faster approach of physical observables to the continuum limit than with the 
unimproved plaquette gauge action. 

In order to satisfy the requirement of chiral symmetry on the lattice, we use the domain- wall formalism 9] for the 
quark action. Adopting the Shamir's formulation jlOl fll| . the action is written as 



S F = -Y;^s{x)D° w {x,y)My), 

xy,st 

D DW = D W + D 5 



(III.2) 
(III.3) 



7m 



2a 



U IXr (x)8{x + fa-y) + 



•7/1 



2a 



Ufa- fi)S(x- fi-y) 
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a 

D 5 st (x,y) = 



M - At 

-8(x-y)6 sU (III A) 



75 jr . 1 + 75 e 

o s +l,t H ^ o s -i 



2d5 ' 2(25 



5(x - y) - — S(x - y)S st , {151.5) 
a 5 



where D w is the ordinary Wilson-Dirac operator in 4 dimensions, M is the domain-wall height which has to be 
adjusted to ensure the existence of chiral modes, e.g., < M < 2 at tree level, and r is the Wilson parameter which 
we choose to be unity. The operator D 5 is the extended part in the fifth direction in which the coordinate is bounded 
byl<s,t<iV 5 . 

Using the chirality projection operators 



P L = —^—> PR = — 7T- > ( IIL6 ) 



quark fields are defined by 



q{x) = P L ^ 1 {x) + P R j} N5 {x), (III.7) 
q(x) = i> N5 {x)P L +i>x{x)P R , (III.8) 

and their mass m/ is introduced as a parameter in the boundary condition in the 5th direction: 

i>N 5 +i(x) = m f aipi(x), ip (x) = m f aip N5 (x). (111.9) 

The operators Qi and Qsub in our numerical simulation are constructed from q and q only, by identifying u,d and s 
with q u , q d and q s . 

Axial vector transformations in five dimensions are defined as 

5ip s (x) = iQ(s)\ a e a s (x)^ s (x), 5$ s (x) = -i^ s (x)Q(s)\ a e a s (x), (111.10) 

where Q(s) = sign(2A^5 — s + 1) and e° (x) is an infinitesimal parameter. This definition leads to the variation 

6q(x) = i l5 \ a e a (x)q(x), (III.ll) 
5q(x) = iq(x)~f 5 \ a e a (x), (11172) 

in terms of quark fields, and the axial-vector current takes the form 
N5 I 

A«(x) = QWg [M*)0- - 7Mx)\ a Mx + £) + i> s (x + A)(i + ^)ul(x)\ a Mx)] ■ (ni.13) 

s=l 

Taking the divergence of A® , one obtains 

V„>4*(x) = 2J 5 Q g (x) + 2m/aP a (III74) 

with 

J 5 a 9 (x) = $N t /i{x)P L \ a il)N t /w{x) - ^ 5 / 2+1 (i)P ;i A''^ 5/2 (x) (III75) 

and 

P a = q(x)X a l5 q(x). (III76) 

The axial vector current A® does not conserve automatically even in the chiral limit m/ — > due to the first term 
Js g on the right hand side. Effects of this breaking term, however, is expected to vanish as N5 — > 00. In practice it 
is necessary to determine the value of N$ for a given set of lattice parameters and a type of gluon action, so that the 
chiral brea king effect due to this term is acceptably small. 

In Refs. |2fl l27j. the chiral property of domain- wall fermion was investigated in detail in the quenched numerical 
simulation. Defining an anomalous quark mass by |26| 

_ (o|Ex^(x^)P 6 (Q)|o) , 
m5qa = (CWlf' ( } 
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the axial Ward-Takahashi identity l|III.14|) yields 

V M ^A«(x)P fe (0)\ = 2a(m f +m 5q )(^2P a (x)P b (0)y (111.18) 

In FIG.n we quote results of m$ q as a function of N5 from Ref. [25I |26| . In the right panel data from the standard 
plaquette gluon action for a -1 w 1 GeV (circles, (3 = 5.65) and a" 1 2 GeV (squares, (3 = 6.0) are summarized with 
two types of exponential fits. The counterparts from the RG-improved gluon action are found in the left panel, where 
(3 = 2.2 and 2.6 correspond to a -1 ps 1 and 2 GeV, respectively. The anomalous quark mass for the RG-improved 
action is an order of magnitude smaller than that for the plaquette action for both a -1 « 1 and 2 GeV. This clearly 
demonstrates the advantage of the use of RG-improved gluon action, which we therefore adopt in our work. 



B. SIMULATION PARAMETERS 



Our numerical simulations are carried out in the quenched approximation at the inverse gauge coupling of (3 = 2.6. 
From the string tension ^fa = 440 MeV j^l, H^, ^3 > this value of (3 corresponds to 

l/a = 1.94(7) GeV, (IH.19) 

which we adopt in our analyses. If we use other quantities such as the rho meson mass or the pion decay constant to 
determine the scale, the lattice spacing is different from the above value, due to the quenched ambiguity as well as 
the scaling violation. We do not include such an ambiguity of a in the systematic uncertainty of our results. 

Denoting the five-dimensional lattice size as iV 3 x N t x N5, we choose the fifth dimensional length to be N$ = 16 
and the domain wall height of the quark action to be M = 1.8. For these parameter choices the anomalous quark 
mass at (3 = 2.6 is given by rns ? = 0.283(42) MeV [26|. We expect this magnitude to be sufficiently small for viability 
of the xPT reduction formulae. Chiral properties of matrix elements will be discussed in detail in Sec. IIV Al 

To investigate the effect of finite spatial volume V — iV 3 , two sizes of lattice given by N s = 16 and 24 are examined, 
in both cases using the temporal size N t — 32. 

We work with degenerate quark masses for u,d and s quarks, and denote the common bare quark mass as rrif = 
m u — met — m s . Matrix elements are evaluated for the bare quark masses mja = 0.02,0.03,0.04,0.05 and 0.06. 
Masses and decay constants of the pseudo-scalar meson calculated on the lattice, which are common for pion and 
kaon, are denoted as tum and /a/ ■ 

Gauge configurations are generated by combining one sweep of the 5-hit pseudo heat bath algorithm and four 
overrelaxation sweeps, which we call an iteration. We skip 200 iterations between configurations for measurements. 
In TABLE [J the numbers of configurations used in our analyses are given. 

We emphasize that we generate gauge configurations independently for each value of raja. This is practically feasible 
since most of the computer time in our runs is spent in calculating quark propagators. A clear advantage is a removal 
of correlations between data at different values of m/ , and hence a more reliable control of the chiral extrapolation 
as a function of my or meson mass squared m 2 M on the basis of \ 2 fitting of data. For error analyses at each mf a 
single elimination jackknife estimation is employed throughout the present work. 

TABLE im shows m 2 M for both sizes of 16 3 x 32 and 24 3 x 32. The intercepts in my and m 2 M are obtained by taking 
a linear extrapolation. Values of m/ in the limit of m 2 M — » are 0.95(62) MeV and 1.09(31) MeV on 16 3 x 32 and 
24 3 x 32 lattices respectively. These values are larger than the value m 5q = 0.283(42) MeV at mj = 0. As pointed 
out in Ref. |26| . the discrepancy between the direct measurement of 7715^ and the estimate from the pion mass is 
largely explained by finite spatial size effects on the pion mass. We use m 2 M as a variable in our chiral extrapolation 
throughout this paper. We have checked that our results remain identical within estimated statistical errors if m/ is 
used in chiral fits. 



C. CALCULATION OF MATRIX ELEMENTS 



In FIG.|2we display the quark line diagrams of three- and two-point correlation functions needed for our simulation. 
Filled squares represent the weak operator Q^p or Q su b located at the site (x, t). Crosses are meson operators. We 
fix gauge configurations to the Coulomb gauge. A wall source for pion is placed at t = and that for kaon at 
t = T = Nt — 1. Quark propagators are solved by the conjugate gradient algorithm, imposing the Dirichlct boundary 
condition in time and the periodic boundary condition in space. The stopping condition is given by 



||( J D + m)-x-6|| 2 <10- 9 ||6|| 2 , 



(111.20) 
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where b is the source vector, x is the solution vector and D is the lattice fermion operator. With this stopping 
condition a precision of better than 0.1% is achieved for arbitrary elements of three-point correlation functions. 

The three-point correlation functions for K — > it matrix elements have the contractions of FIG. 0a), (b) and 
(d). For calculating the I — amplitudes /ir + K + \ both the figure-eight contraction of (a) and the eye 



contraction of (b) are needed, while for the 1 = 2 amplitudes (^n 
Writing 0(t) = y X) x ^( x j *)j we extract the matrix elements from calculation of the ratio of form 



{" 


v+(T)Q { !\t){K+y{Q) 




(0\n+(T)Mt)\0) (0\A i (t)(K+ 


)t(0)|0) 



T>t>l 



K + \ only the figure-eight contributes. 



(111.21) 
(III. 22) 



K+ 



(tt+ |A 4 |0> (0\A 4 \K^ 



IC 



We note that a local current A^(x) = q(x)j^-f5q(x) is employed in the denominator rather than the conserved current 
given in (|III.13(I in order to match with the local form of the four-quark operator in the numerator. 

The contractions in FIG.[3fc) show the K° — > (vacuum) annihilation matrix elements from which the parameters 
OH in the %PT reduction formulae (|II.50J| are obtained. If d and s quarks arc non-degenerate, these parameters are 
easily obtained from the ratio of propagators: 



(111.23) 



In the limit of degenerate quark masses, which applies to our numerical simulation, some care is needed. From the 







J t— >oo 


<° 


Qi 0) 




(0 


Q S ub(t)(K°)H0) 


o> 


(0 


Qsub 





definition of Q su b <|II.41(1 and the fact that CPS symmetry gives (0 \Qi\ K } 



0, we derive 



= — lim 

m s — >md 



(» 






m s > in d 


(m s 


- m d ) 


(0\s~/ 5 d\K°) 




K° 


> 














md 



(0|s 75 d|K ) 



(111.24) 



(111.25) 



The derivative acts both on the operator and on the kaon, and hence there are two contributions as shown in 

tgator is obtained t 

Y J G{x,z)G{z,y). (HI.26) 



FIG. EIc). The necessary derivative of the quark propagator is obtained through 

dG(x,y) 



dm 

z 

To calculate the quark loops that appear in the eye and annihilation contractions, we employ the random U(l) noise 
method. We generate £w (x) = e l6tyX ^ (j = 1, • • ■ , TV) from a uniform random number 9(x) in the interval < 8 < 2n. 
In the limit N — > oo, we have 



i N 



Therefore, calculating quark propagators with (x) as the source, 



(x) = Y,(D + m)-\x,x')^{x% 



(111.27) 



(111.28) 



we find 



x) — > (D 



~ l {x,x) 



(111.29) 



as the quark loop amplitude for each gauge configuration. 
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In our calculation, we generate two noises for each spinor and color degree of freedom, i.e., 2 x (#color) x (^spinor) = 

24 noises for each configuration. In FIGs. |3| and0] we show propagator ratios for the an d Qq operators, and 
those for oti and a§. The horizontal lines indicate the values extracted from a constant fit over t — 10 — 21 and the 
one standard deviation error band. Here correlations between different time slices are not taken into account for the 
fit. Instead errors are estimated by the jackknife method. We observe reasonable signals, which show that 24 noises 
for each configuration we employ is sufficient to evaluate the quark loop amplitude. From (|III.22() . the xPT reduction 
formulae derived in Sec. Ill Cl and III Dl are converted to the following forms at the lowest order of xPT: 



for z = 1, ... ,6,9, 10 



i( 2 ) 



K 



V2fA m K ~ m l) x 
V2f w {m 2 K -ml)x 





Qf 


sub 


K+\ 


(tt+ 


\A 4 1 


0) (0L4 4 |- 


K+) 



;vr+|A 4 |0) (0|Ai|A-+)" 



for i = 7, 8 (I = 0,2) 



7T + 7T 



Q 



(I) 



A 



= -V2U 



m M x 



Q 



a) 



(vr+|A 4 |0) (0\JU\K+y 



(111.30) 
(111.31) 

(111.32) 



where we set pk = {iniM, 6) and p^ = (—irriM, 6) for K + — > ti + matrix elements on the right hand side. We identify 
/m with /, and assign to it the physical value of /„■, since fu agrees with f v in the chiral limit. On the other hand, the 
meson masses m 2 K and m\ in l|III.30|) and (|III.31|I represent the experimental values since they arise from the physical 
K — > 7T7T matrix elements. All of the experimental values used in our calculation are summarized in appendix [5] We 
emphasize that these formulae are valid to the lowest order in xPT. If higher order corrections are small, the right 
hand sides of eqs. I|III.30I) (|III.32|) should depend only weakly on the lattice meson mass m\ 
The two-pion states in the isospin basis are decomposed as 

I Mo) = sj\ K+Tr^ + yllTrV), 
Therefore, matrix elements in this basis are given by (7r + 7r~ K°\ times constants : 



(111.33) 
(111.34) 



)(2) 



((^olQ.I^ ) = 

((7T7r)2 \Qi\ K°) = Vs^ir+n 
We use a shorthand notation 

(Q. i ) / = ((^)/|g 4 |A°), 7 = 0,2 
for the matrix elements in the isospin basis hereafter. 



K° 
K°), 



(111.35) 
(111.36) 

(111.37) 



D. SUBTRACTIONS IN AI = 1/2 MATRIX ELEMENTS 



According to 1III.30|) the contribution of the unphysical operator Q su b has to be subtracted for calculating 



the AI = 1/2 matrix elements. FIG. El shows the original matrix element 



Q 



(0) 



K + ) (circles), the sub- 



m 



K 



traction term — (n + |Q su b| K + ) (diamonds), and their sum (squares), multiplied with a factor V2f 7I ( 
m 2 )/ (n + \A4 1 0) (0 |^4 1 K + ) for conversion to the K — > nn matrix elements (see HIII.30|) '). The left and right columns 
correspond to the spatial sizes 16 3 and 24 3 respectively, and the upper and lower rows exhibit the data for and 
Qg ^ as typical examples. These matrix elements play a dominant role in the AI —1/2 rule and e'/e as we see in later 
sections. The numerical details of subtractions for all of the relevant operators Q f for i — 1,2, 3, 5, 6 are collected in 

TABLE nm 
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We observe that the subtraction term represents a crucial contribution in the physical matrix element. In the case 
of the subtraction term is twice larger than the original matrix element and opposite in sign. Thus, the physical 
matrix element is similar in magnitude but flipped in sign compared to the original matrix element. 

For the case of the subtraction term almost cancels the original matrix element so that the physical matrix 
element is an order of magnitude reduced in size. Nonetheless, as one can see from inspection of TABLE 11111 the 
physical matrix elements are well determined with errors of 10-20%. 

These results show that the subtraction plays a crucial role in calculations with the reduction method. Numerically 
this procedure is well controlled in our case. 



E. RENORMALIZATION AND RG-RUNNING 



Throughout this paper, the renormalization of the operators and the RG-running of the matrix elements are carried 
out within the perturbation theory in MS scheme with NDR. 

The physical K — > ttt: amplitudes in the isospin basis Aj are given by 

r io 

Aj = -%VM Wiifi) (Q l )f S (^), (111.38) 

v i— 1 

where we set Si = since our calculation at the tree level of xPT does not incorporate the effect of the final state 
interaction; this effect begins from the next to leading order of xPT. The Wilson coefficient functions have a form 

Wi(jj) = Ziiri+T-Viifj,) (111.39) 

where yi are non-vanishing only for i = 3, • • • , 10 and r = — (V t * s Vtd) / {V* s V u a) is a complex constant. With our choice 
of scale fi — m c — 1.3 GeV, the functions Zi(m c ) are negligibly small for i = 3, • • • , 10 |4lj | . 

The coefficient functions yi(n) and at m c = 1.3 GeV have been calculated for several values of the QCD 

parameter A— . We employ A— = 325 MeV for our main results, and also consider A— = 215 MeV and 435 MeV 

^ _ MS MS _ ' _ _ MS 

to examine the magnitude of systematic error. The choice of the central value is motivated by recent phenomcnological 
compilations of the strong coupling constant, e.g., Ref. 42] quotes A^ = 296^44 MeV corresponding to a^ s (M z o) — 

0.1184(31). We list the values of coefficient functions we use in TABLE IIVI The experimental parameters are 
summarized in Appendix iBl 

To calculate the renormalized matrix elements in the MS scheme (Qi)r (jii), we first translate the lattice values 
into the renormalized ones at a matching scale q*: 

(Qi)T(q*) = Zijtfa) (Qj)^ tt (l/o). (111.40) 

This step is carried out using the renormalization factor calculated to one- loop order of perturbation theory |43l 
l44l l45l l4t|. The detailed form of the one- loop terms and explicit numerical values for q* = 1/a in quenched QCD, 
appropriate for our case, are given in Appendix IO 

The next step is to evolve the renormalized matrix elements from the scale q* = 1/a to /1 = m c using the 
renormalization group, and combine them with the Wilson coefficient functions Wj(/x). The RG-evolution of the 

matrix elements (Qi)^ S (/i) is inverse to that of the coefficient functions Wi(fx), i.e., 

Wi(jn) = Uin^^ijWjOn), (111.41) 

(Qi)P(An) = [U-^utoftiiiQijfte). (111.42) 

Perturbative calculations of U(m c ,q*) at the next-to- leading order are available [4l|. In Appendix IU1 we adapt the 
known results to calculate the numerical values of the evolution matrix for our case in which /ii = m c = 1.3 GeV 
and fi2 = 1/a = 1.94 GeV. The evolution may be made either for quenched QCD or for Nf — 3 flavors corresponding 
to u, d and s quarks, depending on the view if the matching at fj, = 1/a is made to the quenched theory or to the 
Nf — 3 theory in the continuum space-time. This is an uncertainty inherent in quenched lattice QCD, and we choose 
the Nt = 3 evolution in our calculation. We have also tested the evolution with quenched QCD, and found that the 
results for hadronic matrix elements do not change beyond a 10-20% level. 

For the coupling constant in our Nf = 3 evolution, we employ the two-loop form 



a hiln-r^- 
Pi ln#- 



(111.43) 
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with A— = 372 MeV, which corresponds to A— = 325 MeV. In order to check systematic errors associated with this 

MS ' ^ MS J 



choice, we also make calculations for A-% = 259 MeV (A^L = 215 MeV) and A— = 478 MeV (A£> = 435 MeV). 

' MS v MS 'MS v MS ' 



IV. RESULTS OF HADRONIC MATRIX ELEMENTS 
A. CHIRAL PROPERTIES OF K -> tt MATRIX ELEMENTS 



As we mentioned in Sec. IIII Bl the RG-improved gauge action provides the advantage that the measure of residual 
chiral symmetry breaking m§ q due to finite N$ is small at a -1 ~ 2 GeV. It is nonetheless desirable to check the size 
of chiral symmetry breaking effect directly for the K — > 7r matrix elements. 



Explicit chiral symmetry breaking, if present, causes mixing of the 1 = four-quark operators with the lower 
dimensional operator sd without quark mass suppression, so that K — > -k matrix elements at = m s = rrif behave 
as 



Q 



(0) 



; sub 



„(0) 



) + 



Pi 



m f )(7r+\sd\K+) + 0(m 4 M ) (IV. 1) 



for (8l,1/j) operators, and 



,(0) 



r 2 



(n+\sd\K+)+0{m 2 M ) 



(IV.2) 



for (8i,8fl) operators. Here pi,~ft and Si are dimensionless quantities which represent magnitudes of residual chiral 
symmetry breaking, and hence are proportional to e~ cNa with some constant c. The matrix element (tt + \sd\ K + ) stays 
non-zero in the chiral limit. Motivated by eqs. (|II.41|I and l|III.14(l . one may consider modifications of the subtraction 
operator such as 



Jb — > (m s + m d + 2m 5q )sd - (m s - m d )s"/ 5 d. 



(IV.3) 



Such modifications, however, will not ensure the complete removal of residual chiral symmetry breaking from the 
matrix elements. 

(2) . . _ 

The 1 = 2 operators Q\ 2 do not mix with the sd operator. Their matrix elements can have constant terms in 

(2) 

the chiral limit, however, due to mixings with dimension 6 operators such as Q 7 g in the presence of chiral symmetry 
breaking. Hence we also consider the chiral behavior of these matrix elements. 

Of the ten operators Qi, we recall that (34,9,10 are dependent operators as shown in Hll.llllTT.13l) . Furthermore, 



3,5,6 



there is an identity Q\ = Q 2 which follows from lII.22lHT.23p . and the 7 = 2 component is absent in the Q 
operators. Thus we only need to examine the matrix elements of Q\ \ 3 5 6 an< ^ Q?'' '• 

FIG. ©shows these matrix elements as functions of m 2 M (GeV 2 ) for the two spatial volumes V = 16 3 (left column) 
and V = 24 3 (right column), adopting the normalization defined by 









(7T+ \M 


0) (0 


A 4 \K+) 



x m M a 



2/ 



X 



{I) 



K+ 



(IV.4) 



For the 1 = channel, three data sets are plotted, corresponding to the original matrix element X^ 1 ' = Q\ L> (circles), 



the subtraction term — o^Qsub (diamonds) and the subtracted matrix element Q 

channel, subtractions are absent and hence X^ = Q^K 
The denominator of I|1V.4|) behaves as 



(-0 



<tt+ |A 4 |0><0|A 4 |^+> = 2f M m 



AI ■ 



ctiQsub (squares). For the 1 = 2 



(IV.5) 



irrespective of whether chiral symmetry holds exactly or not. The advantage of our normalization is that the coefficient 
of the m 2 M term of the ratio is directly related to the K° — > ir + ir~ matrix elements. An alternative normalization is 
provided by the ratio 



(IV.6) 





Af> 




(n+\P 


0) (0 


P\K+) 
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where P = 075 g is the pseudo scalar density. This method avoids the use of measured values of pion mass, but it 



loses the straightforward relation to the physical matrix elements. We use the normalization (|IV.4(I in our analyses. 
We have checked, however, that the conclusion remains unchanged even if j|iV.6Jl is employed instead. 
For chiral extrapolation we consider an expansion of the form 

2 

^2- (t + X\ I] K+) = o + ai m 2 M + a 2 {m 2 M f + a 3 (m 2 M ) 2 lnm 2 M + a 4 (m| f ) 3 + • • • . (IV. 7) 



2/ 



M 



Chiral extrapolations using the first three terms are indicated by the solid line in each panel of FIG. |H| The fit 
parameters are summarized in TABLE Ivl The results for the intercept ao in the chiral limit are consistent with zero 
within the fitting errors except for the 1 = 2 operator Q± for the volumes V = 16 3 (1.8er) and 24 3 (4er), the 1 = 



subtracted operator — c^Qsub for V = 16 3 (1.4er) and 24 3 (2.7<r), and the subtraction term for the i = 1 operator 
—ctiQsub for V = 24 3 (2.8ct). Since no systematic tendency that the intercepts become larger for smaller volume is 
observed, it is unlikely that the non-zero intercepts of these matrix elements are caused by the finite spatial size effect. 
Indeed even an opposite tendency that the intercept becomes larger for larger spatial volumes is observed. 

The absence of systematic trend in our data suggests the possibility that non-zero intercepts observed for some of 
the matrix elements are artifacts of the long extrapolation in m 2 M . To test this point, we attempt a fit with a cubic 
polynomial of form a\m 2 M + a2(m 2 l/ ) 2 + a^mj^) and a form with chiral logarithm given by a\m 2 M + a2(m| f ) 2 + 
a 3 (m 2 M ) 2 lnm^j, both having a built-in chiral behavior of vanishing at m 2 M = 0. We show the former fit curves by 
dashed lines in FIG. Eland the fitted parameters in TABLE IvTl Numerical results of the chiral logarithm fit are given 
in TABLE |VIII The fit curves are similar to those of the cubic fit. Both functions provide good fit of data with 
reasonable x 2 /dof. 

Let us try to analyze the chiral behavior of / = matrix elements in terms of mixing with the sd operator as given 
in (|IV.lJl . The existence of the constant /3j can be detected from the chiral limit of the matrix elements. On the 
other hand, separating the contribution of 7$ and Si from the physical ones would require results at different N5. We 
leave such an investigation for future studies, and assume that the latter contributions are negligible. We also ignore 
mixings with the dimension 5 operator sa^F^d since their contributions are subleading in 1/a. 

We estimate Pi from the values of ao obtained in the chiral fit of the matrix elements for the subtracted operator 
Qf^ — a^Qsub given in TABLE For this purpose, we repeat the calculation of (|IV.4|I for = sd, and extract 



\sd\K+)a ^ 2 2 _ 1 



x m M a = 



(ir+\sd\K+) (IV. 



(^+|^4|0)(0|A 4 |^+)a 4 2f 2 M a 

where powers of a are supplied to absorb dimensions of matrix elements. We then fit the results to a quadratic 
polynomial bo + 61 to 2 / + 62 (^m) 2 ■ The numerical values of i|IV.8|) are given in Table lyTTTl and the result for bi are 
given in Table Ivl Normalizing with m^ q = 0.283 MeV to take into account the e~ cN& dependence expected for Pi, 
one has 

A 00 1 (IV.9) 



m 5q a b m 5q a 

In the case of V = 24 3 , the results are (3i/(m 5q a) = 0.9(1.1) for i = 1, -0.91(87) for i = 2, 0.8(4.2) for i = 3, -4.7(4.4) 
for i = 5, and —21.6(8.1) for i = 6. Except for the i = 6 operator for which the coefficient is exceptionally large, we 
find values consistent with zero within the errors. 

The analyses described here do not show strong evidence for the effect of residual chiral symmetry breaking in the 
K — > 7r matrix elements. Although more data at smaller quark masses will be needed for the definite conclusion, 
we conclude here that our results for the matrix elements is consistent with the expected chiral behavior within the 
statistical precision of our data. Therefore, for the chiral extrapolation in the rest of this article, we employ the 
cubic polynomial without constant term for the central value and use the form with chiral logarithm to estimate the 
systematic uncertainty. Since non-zero intercepts beyond statistical errors cannot be excluded for some of the matrix 
elements, we examine possible effects of the residual chiral symmetry breaking to the physical matrix elements in 
SecEl 

Let us also make a comment on comparison of lattice data with predictions of quenched chiral perturbation theory. 
For the 1 = channel, data for more values of mj are required for such a comparison because of the presence of a 
number of unknown parameters as well as a new term of form b\m 2 M lnm 2 ^ in the predicted matrix elements |47l |. On 
the other hand, quenched chiral logarithm terms are absent for the 1 = 2 matrix elements governed by the (27 l, 1r) 

operator, and the ratio 03/01 for is predicted to be 03/01 = — d> / [IQtt 2 f 2 ) = —2.180 GeV -2 . We observe in 
TABLE IVIII that the fitted value agrees in sign but 3 to 4 times smaller in magnitude than the prediction, e.g., 
a 3 /oi = -0.58(10) GeV~ 2 on a 24 3 x 32 lattice. 
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Quenched chiral perturbation theory makes the same prediction for the coefficient of the logarithm term of the chiral 
expansion of Bx as it is governed by the same operator in xPT. For this case, similar discrepancies of lattice results 
from t he p rediction is found for the case of the staggered fermion action [2jj as well as for the domain wall fermion 
action 28] . A possible explanation for these large discrepancies is that higher order corrections in (quenched) xPT are 

(2) 

non-negligible at quark masses employed in the current simulation. Indeed we have confirmed that data for Q\ can not 
be fitted by the form aim 2 M + a2(m 2 M ) 2 + az{m 2 M ) 2 \nm 2 M + a4(m 2 M ) 3 with 03/01 = —2.180 GeV~ 2 fixed. The complete 
form in %PT to this order, aim 2 M + a 2 (m 2 f ) 2 + a?,(rn 2 M ) 2 lnm 2 M + a^mj^) 3 + a^(m 2 M ) 3 lnm 2 ^ + a§{m 2 M ) 3 (lnm^) 2 , 

unfortunately, can not be employed for our data calculated only at five values of quark masses. Understanding the 

(2) 

small value of 03/01 for Q\ ' require further studies. 



B. PHYSICAL VALUES OF HADRONIC MATRIX ELEMENTS 

We tabulate the values of all the K -> tttt matrix elements in TABLEs|IX] (for 16 3 x 32) and|X|(for 24 3 x 32). The 
upper half of each table lists the bare lattice values, (Qi)'/*', and the lower half the physical values, (Qi)^ 5 , obtained 

through matching at the scale g* = 1/a followed by an RG-evolution to /1 = m c . Note that (Q^-q)^ become non-zero 
due to the RG-evolution which breaks the isospin symmetry in the presence of the QED interaction. The two sets of 
numbers do not differ beyond a 10-20% level except for ($5,6,7,8)0; ^ or wmcn the difference amounts to 30-40%. The 
latter situation arises from a larger magnitude of mixing of order 5-10% among the Q 5 ° 6 7 8 operators compared to 
the other operators which are typically less than 5%. In the following, the superscript "MS" will be omitted unless 
the confusion may arise. 

In TABLE IXll we illustrate the magnitude of uncertainty due to the choice of q* by comparing the values of physical 
hadronic matrix elements (Qi}j (m c ) for the choices q* — 1/a and q* — n/a at rrif — 0.02 on a 24 3 spatial volume. 
One finds that the difference is at most 20-30%. 

In FIG.[7|wc plot the physical matrix elements for the AI — 1/2 amplitudes (Qi) (i = 1, • • • , 6, 9, 10) as a function 
of m 2 M . These eight matrix elements involve the subtraction of unphysical effects. The open and filled symbols indicate 
the data from V = 16 3 and 24 3 volumes, respectively. Within the statistical errors at each mt and the fluctuation for 
different values of m/, both of which are larger for the smaller spatial size 16 3 , the data from the two spatial volumes 
do not show indications of presence of finite size effects. 

The remaining matrix elements (Q7,s)o f°r the AI = 1/2 amplitude, which do not require the subtraction, are 
shown in FIG. [SJ These matrix elements are well determined and exhibit a clear m 2 M dependences. 

The matrix elements for the AI = 3/2 channel given by (Qi) 2 = {Qi)% and (Qi ! &) 2 are plotted in FIG. El Their 
statistical quality and m 2 M dependence are similar to those for {Qt^q. 

As discussed in Sec. II V AI for extracting the values in the chiral limit, we adopt a quadratic polynomial form, 

(Qih = + + (iv.io) 

In addition we also employ the chiral logarithm form, 

(Qi)i = Zo+tim 2 M + &m 2 M lnm 2 M . (IV.ll) 

In TABLEs IXIII and IXIIII results from these chiral extrapolations are summarized with the values of % 2 /dof. The 
differences between two types of fits should be taken as a measure of systematic error. For (Qe)o, one observes in 
FIG. [7| an exceptional behavior of the data at m/ = 0.02. An additional chiral extrapolation excluding this quark 
mass is hence also made for comparison and the fit lines indicated in the figures are obtained. 

C. B parameters 

1.3 GeV into B parameters defined by [4l| 

(IV.12) 
(IV. 13) 
(IV. 14) 



We convert renormalized hadronic matrix elements at /1 — m c = 



b{ 1/2) 


= -|(Qt: 


R (l/2) 




B (l/2) 


= !<Q 3 ) 



16 



B, 



B, 



B, 



D 



(1/2) 
5 


3 (O 


o . 


(1/2) 


= y {q& 




(1/2) 




(Qr) 






(1/2) 




(Qs) 



(3/2) 



B, 



(3/2) 



D 



(3/2) 



(Qi 



4V2X ' 
(Qr) 2 

W2 Y+ 72 X 



(IV.15) 
(IV.16) 
(IV. 17) 

(IV. 18) 
(IV.19) 
(IV.20) 

(IV.21) 



where 



/if — fit 



X = V3U(m 2 K -ml) ,Y = -4^3 



'■K 



fx 
K 



(IV.22) 



We summarize the values of B parameters in the chiral limit obtained by the fit with quadratic polynomial or chiral 
logarithm in TABLE IXTVl Quark masses and other parameters used in the calculations are given in Appendix iBl 

Let us compare our values of B parameters with typical ones quoted in phenomenology (see e.g., [41|). For the 
B parameters important for the AI =1/2 rule, the experimental value of KeA 2 indicates B[ N p R (m c ) = 0.453 

■ 0.4 - 0.5 is consistent. On the other hand, our results 



with Alii = 325 MeV, with which our value B {3/2) (m 

MS ' 1 v 



B{ 1/2 \m c ) 



9 and B 



(1/2) 



(m c ) « 3 - 4 are smaller than B 

(1/2) 



(1/2) 



(m c ) 



15 and B 



(1/2) 
2.NDR 



m c ) = 6.6 needed to 



(1/2), 



explain the experimental value of ReAo. For the parameter B\ 1 ' relevant for the direct CP violation, the largest of 
our estimate B, 

fig 1 / 2 -* = 1 in the l/N c approach, while B, 



(m c ) w 0.3 from the 4-point fit of the data from the 24 3 spatial volume is still much smaller than 
§ (m c ) ~ 0.9 is comparable to B^ 2 ^ = 1 again in the 1/N C approach. In 
general the B parameters for 1 = are smaller than usual estimates. 



3 (3/2) 



{/j, = 2 GeV, NDR) 



Previous studies gave By 
2 GeV,RI(MOM)) 

( 3 /2) („ _ 2 GeV, NDR) = 0.80(9) |2 



0.58(7) &ndB ( 8 3/2) (n 



0.38(11) and Bf /2) {^i = 2 GeV, RI(MOM)) 



2 GeV, NDR) = 0.81(4) 20], B { 7 3/2) (/x 



and B^ /z, {^ 

? (3/2 



0.77(9) [21|, B { 7 3/2) (fi = 2 GeV, NDR) = 0.58(9) 



from quenched lattice QCD, and B 7 ° /Z >(fi = 2 GeV, NDR) = 0.55(12) 



and B s ' ' (/i = 2 GeV, NDR) 



values are B 



(3/2) 



1.11(28) from dispersive sum rules where m s + md = 100 MeV is used 48]. Our 
^/ 2 \^ = 1.3 GeV, NDR) = 0.92(4) on a 24 3 x 32 lattice in 



1.3 GeV, NDR) = 0.62(3) and B 8 



broad agreement with the above. Note that the scale /i is different between our results and those of other studies. 



V. PHYSICAL RESULTS 
A. AI = 1/2 RULE 

The real part of Aj relevant for the A7 =1/2 rule is written as 



ReAj = ^\V ud \ ■ \V US \ 



10 



^ Zi{m c ) (Qi) j (m c ) + (Re t) ^ Vi{ m c) (Qi)j (m c ) 



i=3 



(VI) 



In TABLE IXVI we list the values of KeAo, KeA 2 and uj 1 = Re^4o/ReA2 for each value of m/ and spatial volume, 
and for the three choices of the A parameter A™ = 325 MeV, 215 MeV and 435 MeV. 



4 for Ag) 

open and filled symbols denote the results from the volume V = 16 3 and 24 3 , respectively. Signals for Re^2 is quite 



FIG. El plots Re^2 (left panel) and ReA (right panel) as functions of m 2 M for Ai-i = 325 MeV. In both panels, 
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clean, while those for Re^4o exhibit more fluctuations. Since both amplitudes show a variation with m 2 M , we need 
to extrapolate them to the chiral limit to extract the physical prediction. Following the analysis in Section II V Al we 
examine two types of fit functions given by 



Chiral extrapolations from the quadratic fit are indicated by solid lines, and those from the chiral logarithm fit by 
dashed lines in FIG. EH 

For the AI = 3/2 amplitude plotted on the left, the extrapolated values show good agreement with the experimental 
value ReA2 — 1.50 • 10~ 8 GcV indicated by the horizontal arrow. On the other hand, the AI — 1/2 amplitude Re^4o 
is small at measured values of quark masses, and only amounts to about 50-60% of the experimental value 33.3 • 10~ 8 
GeV even after the chiral extrapolation. 

A breakdown of the amplitudes into contributions from the ten operators Qi with i — 1, • • • , 10 is illustrated in 
FIG. ^] for nifa — 0.03. The histograms for the V = 16 3 and 24 3 cases are shown by dashed and solid lines, 
respectively. The horizontal lines with statistical errors indicate the total amplitude, the dashed and solid lines 
corresponding to V = 16 3 and 24 3 . An apparent absence of contributions from the operators with i = 3, • • • 10 is due 
to the small value of the parameter Re r w 0.002; the real part of the decay amplitudes is determined by the matrix 
elements {Qi) I and {Q2) j, with the latter providing the dominant part. 

The ratio lo^ 1 = ReAo/ReA2 is shown in FIG. El Reflecting an insufficient enhancement of the AI — 1/2 

amplitude, it only rises to about half of the experimental value uj~ 1 » 22. The situation hardly changes for A^jL = 
215 MeV or 435 MeV, for which the amplitudes shift by about 5-10% (see TABLE |XVj . We collect chiral fit 
parameters for the case of larger spatial volume V = 24 3 in TABLE IXVll 
Altogether we find 



The central values are taken from the result on a 24 3 x 32 lattice from the quadratic polynomial fit with ALA. = 325 
MeV. The first error is statistical, the second one is an estimate of uncertainty of chiral extrapolation using the chiral 
logarithm fit, the third one is finite-size variation estimated by the change of value for for the V = 16 3 lattice, the 
fourth one, associated with renormalization, is estimated as the largest variation under changes of A^-, q* and the 

RG-running. If the chiral symmetry breaking term, £_i/m^, is included in the chiral fit (|V.2fl . a non-zero value of 
£_i beyond the statistical error is obtained only for ReA2, resulting in a 60% increase of the value of Re^- The 
disagreement from experiment becomes worse in this case. The scaling violation and the quenching error, which 
cannot be estimated in our calculation, are not included in our systematic uncertainty. In particular, the physical 
scale of lattice spacing set by the string tension in this paper may differ by about 10-20% from scales determined 
by other physical quantities due to the quenched approximation. This uncertainty is not included in the above error 
estimate. 




(quadratic polynomial) 



(chiral logarithm). 



(V.2) 




(V.3) 
(V.4) 
(V.5) 



B. DIRECT CP VIOLATION (e'/e) 



The formula (|I.4(1 for e'/e can be rewritten as 



e'/e 



lm(V*V td ) pW 2 ) - p( 3 / 2 ) 



(V.6) 




(V.7) 




(V.8) 



where 



(V.9) 



r = 



2|e|Re^ 
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and the parameter £l v +rj' = 0.25(5) reflects the isospin breaking. Since the AI = 1/2 rule is only partially reproduced 
with our data, we employ the experimental values for HeA , w and e as input. 

In FIG.EHour data for P (3/2) (left panel) and P (1/2) (right panel) calculated with = 325 MeV are plotted as 

a function of m 2 M . Results for e'/e are shown in FIG. El Since P^/ 2 ) is smaller than p( 3 / 2 ) in our data, e'/e tends 
to be negative. 

A breakdown of p( 3 / 2 ) and P^ 1 / 2 ) into contributions from the operators Qi (i — 3, • • • , 10) is displayed for the case 
of m/et = 0.03 in FIG. El where dashed and solid lines denote data from V = 16 3 and 24 3 , respectively. This figure 
demonstrates that {Q%) 2 an d (Qs)n are respectively dominant in p( 3 / 2 ) and P^ 1 / 2 ) as usually considered. However, 
the matrix element of (Qe)o i s too small; if the experimental value of e'/e is to be reproduced by a change of this 
matrix element, it has to be increased by about a factor of five. 

Numerical values of pt 1 / 2 ), p( 3 / 2 ) and e'/e for each m/ are summarized in TABLE IXVIII In addition to the 

features of data discussed above, we observe that changing the A parameter from A^ = 325 MeV to 215 MeV 

decreases Pt 1 ^ by 20% and P< 3 / 2 ) by 25%. Employing A^L = 435 MeV leads to an increase by similar percentages 
for the two functions. Therefore, the trend toward a negative value of e'/e is not altered. 

If we make a quadratic chiral extrapolation we find e'/e = —7.7(2.0) x 10~ 4 with % 2 /dof = 1.75 on a 24 3 x 32 
lattice. Including the chiral symmetry breaking term £,-i/m 2 M in the fit changes this value to +30(20) x 10 -4 with 
X 2 /dof = 0.0015. The small x 2 indicates that more data points, in particular data at smaller masses, are necessary 
to constraint the fit parameters well. The existence of large uncertainties associated with the possible presence of 
the chiral breaking term, and also a subtle quenching effect mentioned below, make it difficult to draw a conclusive 
estimate of e'/e. 

Recently, Golterman and Pallante pointed out that the relation between K — > tt and K — ► ttit matrix elements in 
chiral perturbation theory should be modified in the quenched theory ^| . We have applied the modified relation to the 
Q^l matrix elements, and found that the effect is large, ranging between 20% to 100% in magnitude. For example, the 
renormalized (Q 6 ) on a 24 3 x 32 lattice increases in magnitude to -0.154(17), -0.182(16), -0.144(11), -0.1238(90), 
—0.0969(72) at mj = 0.02, 0.03, 0.04, 0.05, 0.06, respectively. (This modification has been tested also in the case 
of the staggered fermion |5(|, and an increase of (Qg)q of a similar magnitude has been observed.) In terms of 

e'/e, the modified relation leads to -1.70(53), -0.53(51), -1.48(32), -2.09(26), -2.85(19) for A__ = 325 MeV. The 
modification increases the value of e'/e, but it is still negative. A complete analysis still remains to be made both in 
the theoretical analyses of the relation in quenched chiral perturbation theory and in numerical simulations. 



VI. CONCLUSIONS 



In this article we have presented results of our investigation into the reduction method in the framework of chiral 
perturbation theory at the lowest order to calculate the K — ► ttit decay amplitudes. The K — ► tt and K — > hadronic 
matrix elements of four-quark operators were calculated in a quenched numerical simulation using domain- wall fermion 
action for quarks and an RG-improved gauge action for gluons to satisfy the requirements of chiral symmetry on the 
lattice. We have seen that the calculation of quark loop contractions which appear in Penguin diagrams by the random 
noise method works successfully. As a result the A J =1/2 amplitudes which require subtractions with the quark 
loop contractions were obtained with a statistical accuracy of about 10%. We have investigated the chiral properties 
required for the K — > 7r matrix elements. If we leave aside we have found no strong sign for the existence of the 
chiral symmetry breaking effect within the statistical precision of our data in the range of quark masses employed in 
our simulations. However, Qg appears to show an exceptionally large chiral symmetry breaking effect compared to 
other channels. It is not clear to us if this is an effect beyond statistical fluctuation. For the definite conclusion on 
this point, more data, particularly at smaller quark masses, will be needed. Matching the lattice matrix elements to 
those in the continuum at /x = 1/ct with the perturbative renormalization factor to one loop order, and running to the 
scale fi = m c = 1.3 GeV with the renormalization group, we obtained all the matrix elements needed for the decay 
amplitudes. Unfortunately the physical amplitudes thus calculated show unsatisfactory features. 

One of the pathologies of our results is a poor enhancement of the AI = 1/2 decay amplitude; the value of KeA^ 
is about 50-60% of the experimental one in contrast to HeA2 which reaches the expected value in the chiral limit. 
Another deficiency is a small value of the AI = 1/2 contribution to e'/e; if we assume that the AI = 3/2 contribution 
has a correct order of magnitude, the AI = 1/2 contribution is too small by about a factor of five to explain the 
experimental value ~ 2 • 10~ 3 . 

The hadronic matrix elements for AI = 1/2 involve significant subtractions. For some of the matrix elements, 
this results in flips of sign and a reduction in the magnitude. Hence insufficient choices of lattice parameters in 
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simulations may lead to sizable systematic errors in these matrix elements. Possible origins of the errors are (i) 
finite fifth dimensional size N 5 of the domain wall fcrmion, (ii) finite spatial size N s , (iii) finite lattice spacing a, 
(iv) quenching effects, and (v) the neglect of the charm quark. Our use of (vi) renormalization factors in one-loop 
order of perturbation theory is another source of error in the renormalizcd matrix elements. Finally (vii) higher 
order corrections in chiral perturbation theory is also a possible source of error. It may well be that the origin of the 
deficiency resides in physical phenomena such as the effect of a resonance which are difficult to take into account once 
the reduction to K — > ir matrix elements is made. 
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APPENDIX A: DECOMPOSITION OF Q 4 'S INTO AJ = 1/2 AND AI = 3/2 PARTS 

Four-quark operators which transform under the irreducible representations of SU(3)l *S> SU(3) r chiral group and 
having definite isospin / = or 2 are given by 



y(2) 
^27,1 


= {sd) L [{uu) L 


- (dd) L ] + 


(su) L (ud) Ll 


(A.l) 


y(0) 
^27,1 


= {sd) L [{uu) L 


+ 2{dd) L - 


' 3(ss) L ] + (su) L (ud) L , 


(A.2) 


^8,1 


= (sd) L (uu) L - 


■ (su) L (ud) 


L, 


(A.3) 


•G(O) 
^8,1 


= {sd) L [{uu) L 


+ 2{dd) L 4 


- 2(ss) L ] + (su) L (ud) L , 
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y(0) 


= {sd) L [{uu) R 


+ (dd) R + 


(ss)n\ , ^3 c , 


(A.5) 


v (o) 

^8.8 


= {sd)L [(uu)r 


- (ss) R ] - 


(su) L (ud) R , y£l c , 


(A.6) 


^8,8 


= {sd) L [(uu) R 


- (dd) R ] 4 


- {su) L (ud) R , y£> B c , 


(A.7) 



where we use the notation of X's and for the Lorentz structure L (g> L and L <g> R. The subscripts "i, j" stand for 
the representation {iliJr) of the operator and the superscript (0) or (2) denotes the isospin. A shorthand notation, 
e.g., (sd)i = S7/i(l — 75 M is employed as in eqs. (|II.17|) - (III.21(1 . and y^ 1 ) c equals y^O with its color summation 
changed to cross the two currents. In terms of these operators the independent local operators are rewritten as 
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(A.8) 

(A.9) 

(A.10) 
(All) 
(A.12) 

(A.13) 
(A.14) 

Therefore, the decomposition of the local operators into A J =1/2 and A J = 3/2 parts is summarized as follows: 
AI = 1/2: 

1 _ _ _ - 

-{s a d b ) L (u b u a )L + 2(s a u b ) L (u b da)L + (s a d b ) L {d b d a ) L ], 



v (o) 

^8,8 



Q? 
Qf 
Qt 



3 

-[-{sd) L {uu) L + 2{su) L {ud) L + (sd) L {dd) L ], 

(sd) L [(uu) L + (dd) L + {ss) L }, 
{s a d b ) L [(u b u a ) L 4- (d a d b ) L + (s b s a ) L ], 
{sd) L [(uu) R + (dd) R + (ss) R ], 
{s a d b )L[(u b u a ) R + {d a d b ) R + {s h s a ) R ], 
^[(sd) L (uu) R - (su) L (ud) R - {sd) L (ss) R }, 

z\\_{s a d b )L{u b u a ) R - {s a u b ) L (u b d a ) R - {s a d b ) L (s b s a ) R }, 

-[(sd) L (uu) L - (su) L (ud) L - (sd) L (ss) L ], 

-j\.(sad b ) L {u b u a ) L - (s a u b ) L (u b d a ) L - (s a d b ) L (s b s a ) L ], 



(A.15) 

(A.16) 

(A.17) 
(A.18) 
(A.19) 
(A.20) 

(A.21) 
(A.22) 
(A.23) 
(A.24) 
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A/ = 3/2: 



Qi 2) = Q { 2 ] = \[{sd) L (uu) L + (su) L (ud) L - (sd) L (dd) L ], (A.25) 

QW = Q« = QW = = o, (A.26) 

Qt 2) = ^[(sd) L (uu) R + (su) L (ud) R -(sd) L (dd) R ], (A.27) 

is = ^[{sadb)L(ubU a ) R + {s a u b ) L {u b d a ) R ~ (s a d b ) L (dbd a ) R ], (A. 28) 



,(2) _ n (2) _ 3 ( 2 ) 



Qtf = 2 Q i ' (A - 29) 



where color indices are understood within each current in the operators with two color traces. The equivalence between 
and Qj is valid due to Fierz rearrangement, hence = follows. 



APPENDIX B: EXPERIMENTAL INPUT PARAMETERS 

We collect the input parameters which were used in our numerical calculation |5ll l52| . 



Quark mass : m u = 5 MeV , = 8 MeV, (B-1) 

m s = 120 MeV , m c = 1.3 GeV, (B.2) 

m b = 4.2 GeV , m t = 170 GeV (B.3) 

Meson mass : = 139.6 MeV , m K = 497.7 MeV (B.4) 

Decay constant : f„ = 92.4 MeV , f K = 113.1 MeV (B.5) 

Coupling constant : a = e 2 /(47r) = 1/129 (at fj, = mw), (B-6) 

G F = = 1.166 • 10~ 5 GeV" 2 (B.7) 
8m^ 

(m w = 80.2 GeV) 

Quantities relevant 

to Kaon decays : ReA = 33.3 • 10 -8 GeV (B.8) 

ReA 2 = 1.50 ■ 10" 8 GeV (B.9) 

\u\ = 0.045 (B.10) 

n, )+ri , = 0.25 (B.ll) 

|e| = 2.280 • 10~ 3 (B.12) 

CKM elements : |Ks| = 0.22 , |Kd| = 0.974 (B.13) 

Im(V*V td ) = 1.3 • 10- 4 (B.14) 

Rer = -Re ( Yzs^L) = 0.002. (B.15) 



APPENDIX C: RENORMALIZATION FACTORS AND RG-EVOLUTION MATRIX 

In this appendix, we summarize the renormalization factors and the RG-evolution matrix, and calculate their 
numerical values for our choice of parameters. Throughout this paper, we employ the perturbative calculation in MS 
scheme with NDR. 

The renormalization formula has the form, 

(Qi) m iq*) = Zfjtfa) (Qf u ) (1/a) + Zf^tfa) (O (1/a), (C.l) 

where, 

Q^=Q 4 + Q 6 - ( ° 3 + ° 5) i (N c = 3: #color) (C.2) 
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is the sum of contributions from penguin operators. Since our matrix elements are obtained in the form of propagator 
ratios, Z 9 and Z pen are also ratios of the renormalization factors Zf. and Zf cn calculated from corresponding vertex 



functions and that of the local axial current Za [43 



The diagonal parts Zf { are given by 



Z 9 - 



7 2 



(C.3) 
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16tt 2 
9 2 



N, 
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z + + Z- 



16tt 2 



16tt 2 



-— ln(q*a) 2 + z\ - v 2 \ 



m 2 c 



■ ln(q*a) 2 + z 2 + v 2 \ 



while for off-diagonal parts, one has 



7 9 — 



16tt 2 



16tt 2 



3ln(q*a) 2 + 



Z 2 - Z\ + V 2 \ - V\2 



167T' 



■N c v 21: 







1,2,3,4,9,10, 
i = 5, 7, 
i = 6,8, 



(t,j) = (1,2),(2,1), (3,4), 
(4, 3), (9, 10), (10, 9), 
z = (5,6),(7,8), 

i = (6,5),(8,7), 
others 



(C.4) 



(C.5) 



Similarly the contributions from the penguin operators [4jj are given by 

g 2 a 



7 pen _ 



167T 2 3 



ln{q* a y + z\ 



(C.6) 



where C 2 = 1, C 3 = 2, C 4 = C* 6 = AT/, C* 8 = Cio = N u — N d /2, C 9 = -1 and C t = for other i with iV>, AT„, iV d being 
the number of flavors, up-like quarks and down-like quarks in Qi's, and zf en are constants. In our calculation, we 
should set Nf — 3, N u = 1, iV<j = 2. Finally the axial vector renormalization constant has the form, 



Z A = 1 



12tt : 



■ZA 



(C.7) 



In the above z± 1 z\ 1 z 2 , v\ 2 , v 2 \ and za are constants depending on the choices of simulation parameters and renor- 
malization scheme. With the use of mean field improvement at one-loop level, we obtain the following values 0] at 
(3 = 2.6 and M = 1.8 for the RG-improved gauge action: 



g 2 ee g^(\/ a) - 2.273 

M = 1.41979 
z A = -4.6930 

z + = -13.612 , z_ = -10.319 

2i = -10.063 , z 2 = -16.125 

V12 = 8 , V21 = 1 

„.,, j 4.494 (fori = 2, 3, 5, 7, 9) 
^ 3.494 (for i = 4,6,8,10) 



(C.8) 

(C.9) 
(CIO) 
(C.ll) 
(C.12) 
(C.13) 

(C.14) 



From the definition of Q pon , Zf cn can be written in the form of a 10 x 10 matrix Z pen , defined as Zf^ n 



X i5 



-zr/N c ,zf 4 -~ i6 

10 x 10 matrix given by 



7-pcn y-pcn 



pen 



, and Zf- = for other j. The renormalization factor can then be summarized as a 



Z g , Z pcn 
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0.9997 -0.0350 

-0.0350 0.9997 -0.0106 0.0318 -0.0106 0.0318 

0.9785 0.0287 -0.0212 0.0636 

-0.0597 1.0739 -0.0247 0.0742 

1.0154 -0.0924 

-0.0247 0.0742 -0.0884 1.0190 















1.0154 -0.0924 

-0.0637 0.9448 

0.0106 -0.0318 0.0106 -0.0318 0.9997 -0.0350 

-0.0350 0.9997 



(C.15) 



For the derivation of the RG-evolution matrix, we start with constructing the renormalization group equation (RGE) 
of Wi(/i)'s, and hence of U(fi, l/a)'s. If we write the renormalization of Qi as = ZijQj where the superscript (0) 
indicates the value at tree level, RGE for Q^s are readily obtained as 



d In /i 



Qi 



1- 



din (i 



(C.16) 



On the other hand, interpreting Wj's as coupling constants in the effective Hamiltonian, renormalization of Wj's is 
possible, W, — Z^Wj, in place of that of Qj's. From the equivalence of these renormalizations, Z c = (Z^ 1 ) T 
follows. Therefore, using (|G.16|I . we obtain 

(C.17) 



W i = 'y? j W j , hence — — E/y(/x, l/a) = (^ T ) ik U k j(fJ,, l/a). 



dhxfj, 



dln/i 



Using the 10 x 10 anomalous dimension matrix 7, defined in (|G.16|) . the RGE for UUiA/a) has been solved for the 
QCD j3 function and anomalous dimension 7 calculated at next to leading order |53l l54j|: 



is{g 2 



Pa 
is{g 2 ) 

The solution at this order is written as 



16tt 2 



01 



(16vr 2 ) 5 



-in 



r(g 2 



(o)«s . (1) /asy 
Ts'-^+'Ysi — ) 



\AttJ 



4tt 



4tt 
(0)1 _ 



Using the matrix V that diagonalize the 7^ , we obtain diagpyj^] = V 1 j(°') T V and G = V 1 7^ 1 ^ T U. Then 

47T 47T 



/'2 



V 



\Ois(pl) 

J = VHV~\ 



Moreover, with = V~ X ^ )T V 



x (0) #1 



G, 



(C.18) 

(C.19) 
(C.20) 
(C.21) 
(C.22) 

(C.23) 

3 

(C.24) 

(C.25) 
(C.26) 
(C.27) 



2tt 



R{Hi,Hz) = —w-V 



K(°)( Ml , M2 ) + i-^^ (1) ( m , M2 ) 



(C.28) 



24 



2/3 M 



(«) 



(o) 



(0) 

as {^2 



2ft 

7$ A*. 



i^l 1) (Ml.M2)ij 



M« = 
^l X) (/ii,M2) = 



(i) 



7_Di 



t (0) 



'as(M2 



\a,s(A*i 

A 76 

as( M i)^ (0) (Mi,M2), 



7™/2/3o 



as(^2) 



7^/2/3o' 



(C.29) 

(C.30) 

(C.31) 

(C.32) 
(C.33) 



where, = fi c = 1.3 GeV, /12 = l/a- 

Using the value of the strong coupling constant af^(i/ a ) = 0.30171 and af § '(1.3 GeV) = 0.39601 with A_- = 372 
MeV, together with 7-functions presented in Ref. [4lj . we obtain the matrix U(m c , l/a, a) given in (ESS^and the 
RG-evolution matrix: 



0.9738 


0.0730 


0.0035 


-0.0003 


-0.0033 


-0.0002 


0.0005 





0.0005 


0.0001 


0.0731 


0.9736 


-0.0024 


0.0149 


-0.0053 


0.0116 


0.0002 





0.0001 


0.0001 








0.9794 


0.1043 


-0.0212 


0.0247 


-0.0002 





-0.0006 


-0.0001 








0.0731 


1.0105 


-0.0186 


0.0306 


-0.0005 





-0.0004 











-0.0083 


-0.0065 


1.0465 


-0.0996 


0.0005 

















-0.0090 


0.0228 


-0.0421 


0.7878 





0.0007 




















0.0002 





1.0349 


-0.0929 


0.0008 




















0.0004 


-0.0367 


0.7602 


0.0001 


-0.0001 








0.0021 


-0.0149 


0.0053 


-0.0116 


0.0009 


0.0001 


0.9750 


0.0731 








-0.0035 


-0.0001 


0.0032 


0.0002 


0.0006 





0.0736 


0.9740 



(C.34) 



In order to check the systematic error associated with the matching procedure above, we also employ an alternative 
procedure in which the RG-evolution is carried out in the quenched theory from /12 = q* to fii = fj, c =1.3 GeV where 
matching to the Nf = 3 theory is made. For the quenched RG-evolution, the two-loop anomalous dimension matrix 

7g is modified according to [54| 



r 

[lS J quenched 



(C.35) 



where = diagfTi, T2, T3, T4, T5] with the 2x2 matrices which are given by 



ri = r 2 = r 5 = 



2N f 

3N C 
2N, 



22N S 

AN f 



2Nj_ 

3 

2N f 
3N C 



22N 



f 



20C F N f 
3 



4Nj_ 



(C.36) 
(C.37) 



Note that Nf = 3 in this case. For the gauge coupling in the quenched theory, we employ a s is (l/a) 
and o^ s (1.3 GeV) = 0.20439 obtained by the 2-loop running with Nf = 0. 
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rrifa 


16 3 x 32 


24 3 x 32 


0.02 


407 


432 


0.03 


406 


200 


0.04 


406 


200 


0.05 


432 


200 


0.06 


435 


200 



TABLE I: Number of gauge configurations, independently generated for each value of m/a, in our numerical simulation. 



16 3 


x 32 


24 3 


x 32 


m/a 


m| f [GeV 2 ] 


m/a 


m|f[GeV 2 ] 


-0.00049(32) 


0.00 


-0.00056(16) 


0.00 


0.00 


0.0059(37) 


0.00 


0.0066(19) 


0.02 


0.2434(26) 


0.02 


0.2445(11) 


0.03 


0.3568(29) 


0.03 


0.3534(17) 


0.04 


0.4741(28) 


0.04 


0.4714(19) 


0.05 


0.5932(29) 


0.05 


0.5957(19) 


0.06 


0.7134(30) 


0.06 


0.7158(20) 



TABLE II: Lattice pseudo scalar meson mass squared m 2 M [GeV 2 ] at each m/a. The x- and y-intercepts are obtained through 
a linear chiral extrapolation. Physical scale of lattice spacing equals 1/a = 1.94 GeV determined by y/a — 440 MeV. 
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16 3 x 32 

m/ a first subtraction total 


24 3 x 32 
first subtraction total 


Qf ] 0.02 -0.0135(44) -0.0134(25) -0.0269(56) 
0.03 -0.0084(29) -0.0133(20) -0.0217(39) 
0.04 -0.0091(21) -0.0107(16) -0.0198(30) 
0.05 -0.0096(16) -0.0085(13) -0.0181(24) 
0.06 -0.0071(14) -0.0073(12) -0.0144(20) 


-0.0028(23) -0.0164(12) -0.0192(28) 
-0.0052(24) -0.0115(14) -0.0167(30) 
-0.0082(14) -0.0092(10) -0.0174(21) 
-0.0056(12) -0.0076(10) -0.0131(18) 
-0.0085(11) -0.00752(82) -0.0160(16) 


Qf 0.02 -0.0410(40) 0.0825(25) 0.0415(46) 
0.03 -0.0419(24) 0.0755(21) 0.0336(29) 
0.04 -0.0392(18) 0.0752(18) 0.0361(22) 
0.05 -0.0375(15) 0.0659(14) 0.0284(16) 
0.06 -0.0346(13) 0.0627(13) 0.0281(14) 


-0.0500(19) 0.0875(13) 0.0375(23) 
-0.0450(19) 0.0823(16) 0.0373(22) 
-0.0434(18) 0.0743(14) 0.0309(16) 
-0.0394(13) 0.0680(12) 0.0286(13) 
-0.0354(11) 0.0636(11) 0.02821(94) 


Qf ] 0.02 -0.130(17) 0.1253(90) -0.005(21) 
0.03 -0.118(10) 0.1107(81) -0.007(14) 
0.04 -0.1137(76) 0.1179(61) 0.004(11) 
0.05 -0.1132(65) 0.1055(50) -0.0077(86) 
0.06 -0.1000(50) 0.1031(47) 0.0032(75) 


-0.1151(81) 0.1256(46) 0.010(10) 
-0.1140(84) 0.1311(53) 0.017(11) 
-0.1198(61) 0.1206(42) 0.0008(80) 
-0.1052(46) 0.1146(43) 0.0094(66) 
-0.1049(44) 0.1051(35) 0.0002(55) 


Qf ] 0.02 1.719(45) -1.743(36) -0.024(24) 
0.03 1.608(37) -1.657(32) -0.048(15) 
0.04 1.591(33) -1.633(30) -0.042(11) 
0.05 1.438(26) -1.482(25) -0.0444(82) 
0.06 1.430(26) -1.465(23) -0.0359(71) 


1.832(24) -1.853(185) -0.022(11) 
1.731(31) -1.768(261) -0.036(11) 
1.593(27) -1.635(249) -0.0420(80) 
1.521(25) -1.553(224) -0.0321(67) 
1.412(23) -1.448(205) -0.0361(53) 


Qf ] 0.02 4.98(13) -5.01(10) -0.025(51) 
0.03 4.66(10) -4.792(91) -0.129(26) 
0.04 4.632(97) -4.721(86) -0.089(19) 
0.05 4.155(78) -4.287(71) -0.132(12) 
0.06 4.121(73) -4.234(67) -0.1129(95) 


5.264(67) -5.350(54) -0.086(22) 
4.960(88) -5.110(76) -0.150(20) 
4.595(80) -4.732(71) -0.137(14) 
4.385(72) -4.496(65) -0.111(11) 
4.087(65) -4.183(60) -0.0957(88) 



TABLE III: Subtraction in K — > ix matrix element ( 7r 



-,(()) 



K + \ for i — 1,2,3,5,6 multiplied with a factor V^/ir(' 



m 2 K 



m%)/ (k + \Aa\ 0^ ^0 | Ail K + ^. The values of the K + — > 7r + matrix element (first), the subtraction term — at (n + |Q S ub| K + } 
(subtraction) and their sum (total) are given in units of GeV 3 . 



MS 


Zl 


Z2 


2/3 


yi 


2/5 


y& 


yr/a 


2/8 /a 


y 9 /a 


2/io jot 


215 MeV 


-0.346 


1.172 


0.023 


-0.048 


0.007 


-0.068 


-0.031 


0.103 


-1.423 


0.451 


325 MeV 


-0.415 


1.216 


0.029 


-0.057 


0.005 


-0.089 


-0.030 


0.136 


-1.479 


0.547 


435 MeV 


-0.490 


1.265 


0.036 


-0.068 


0.001 


-0.118 


-0.029 


0.179 


-1.548 


0.664 



TABLE IV: Wilson coefficient functions 0. 
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16 3 x 32 






24 3 x 32 










ao 


ai[GeV" 2 ] 


a 2 [GeV" 4 ] 


X 2 /dof 


ao 


ai[GeV" 2 ] 


a 2 [GeV~ 4 ] 


X 2 /dof 






-0.007(38) 


-0.09(17) 


0.04(17) 


0.63 


0.017(24) 


-0.10(11) 


0.01(12) 


1.88 




CflQsub 


0.004(25) 


-0.17(12) 


0.15(13) 


0.19 


-0.041(15) 


0.037(74) 


-0.064(81) 


0.06 


q5 0) 


— aiQsub 


-0.007(51) 


-0.24(24) 


0.17(24) 


0.12 


-0.024(31) 


-0.06(15) 


-0.06(16) 


1.16 




Q ( 2 0) 


0.021(34) 


-0.51(15) 


0.25(15) 


0.07 


-0.002(23) 


-0.50(11) 


0.26(12) 


0.35 




CttQsub 


0.000(26) 


0.82(13) 


-0.37(13) 


2.43 


0.017(18) 


0.802(92) 


-0.37(10) 


0.71 




— a 2 Qsub 


0.019(40) 


0.31(18) 


-0.13(18) 


2.09 


0.024(25) 


0.26(11) 


-0.06(12) 


1.40 






0.02(15) 


-1.36(64) 


0.59(65) 


0.42 


0.063(90) 


-1.43(44) 


0.58(47) 


0.92 




O^Qsub 


0.0014(95) 


1.12(46) 


-0.31(49) 


0.89 


-0.089(59) 


1.71(30) 


-0.92(33) 


0.34 


Qf ] 


— a3Q su b 


0.02(19) 


-0.19(86) 


0.23(89) 


0.58 


-0.02(11) 


0.27(56) 


-0.32(60) 


0.87 






0.37(48) 


14.7(2.3) 


-3.7(2.5) 


2.56 


0.27(34) 


16.8(1.7) 


-6.4(1.9) 


0.48 




C*5<3sub 


-0.15(42) - 


-16.1(2.0) 


4.7(2.2) 


2.56 


-0.11(29) 


-17.8(1.5) 


7.0(1.7) 


0.40 




— asQsub 


0.020(20) 


-1.27(90) 


0.95(92) 


0.17 


0.13(12) 


-0.85(57) 


0.51(61) 


0.48 




Qt ] 


0.8(1.4) 


44.1(6.6) 


-12.4(7.1) 


2.96 


0.86(97) 


47.6(5.0) - 


-17.2(5.6) 


0.31 




aeQsub 


-0.1(1.2) - 


-47.8(5.8) 


14.8(6.3) 


2.61 


-0.19(85) - 


-52.1(4.5) 


21.0(5.0) 


0.42 




— QeQsub 


0.053(38) 


-3.0(1.6) 


1.7(1.6) 


2.27 


0.59(22) 


-4.0(1.0) 


3.3(1.1) 


0.62 




Q< 2) 


-0.0023(13) 


0.0727(64) 


0.0178(68) 0.19 


-0.00264(65) 


0.0751(33) 


0.0140(37) 0.28 






bo 


h [GeV" 2 ] 


b 2 [GeV~ 4 ] 


X 2 /dof 


bo 


bi [GeV- 2 ] 


b 2 [GeV" 4 ] 


X 2 /dof 




sd 


-170(11) 


116(46) 


-64(45) 


2.21 


-186.2(4.0) 


151(19) 


-82(19) 


3.72 



TABLE V: Fit parameters for the chiral extrapolation of the K — > 7r matrix elements defined by ljIV.4(l which should vanish in 
the chiral limit. The parameters (ao, ai, a 2 ) are determined by the fit function ao + aimfj + a 2 (m 2 v/ ) 2 . 



16 3 x 32 







ai[GeV" 2 ] 


a 2 [GeV" 4 ] 


a 4 [GeV- 6 ] 


X 2 /dof 


ai[GeV" 2 ] 


a 2 [GeV" 4 ] 


a 4 [GeV- 6 ] 


X 2 /dof 




Qi 0) 


-0.10(11) 


0.01(44) 


0.06(41) 


0.64 


0.004(67) 


-0.18(28) 


0.10(27) 


2.07 




QlQsub 


-0.159(72) 


0.14(30) 


-0.01(29) 


0.20 


-0.271(40) 


0.64(18) 


-0.49(18) 


0.12 




— aiQ su b 


-0.28(15) 


0.21(60) 


0.00(57) 


0.13 


-0.257(87) 


0.42(37) 


-0.36(36) 


0.98 






-0.37(10) 


-0.06(39) 


0.21(37) 


0.10 


-0.499(63) 


0.23(27) 


0.02(28) 


0.35 




a 2 Qsub 


0.826(75) 


-0.39(31) 


0.02(30) 


2.43 


0.942(47) 


-0.72(21) 


0.26(22) 


0.43 


Of 


— Q2Qsub 


0.46(12) 


-0.46(45) 


0.23(42) 


2.07 


0.461(67) 


-0.56(28) 


0.37(27) 


0.98 






-1.09(43) 


-0.2(1.7) 


0.6(1.6) 


0.35 


-1.00(25) 


-0.3(1.0) 


0.6(1.0) 


1.02 




Ct3Qsub 


1.20(27) 


-0.4(1.1) 


0.1(1.1) 


0.90 


1.07(16) 


0.51(71) 


-0.97(71) 


0.57 




— QsQsub 


0.05(55) 


-0.4(2.2) 


0.5(2.1) 


0.55 


0.10(32) 


0.1(1.3) 


-0.3(1.3) 


0.87 




Qt ] 


17.7(1.4) - 


-11.1(5.7) 


5.6(5.6) 


2.36 


18.82(88) 


-11.0(4.0) 


3.3(4.1) 


0.50 




C*5<2sub 


-17.5(1.2) 


8.5(5.0) 


-3.0(5.0) 


2.44 


-18.68(76) 


9.2(3.5) 


-1.6(3.7) 


0.38 


qT 


— QsQsub 


0.08(59) 


-1.9(2.3) 


1.9(2.2) 


0.23 


0.03(34) 


-1.3(1.4) 


1.2(1.4) 


0.68 






50.9(3.9) - 


-29(16) 


13(16) 


2.80 


54.0(2.5) 


-32(12) 


10(12) 


0.34 




Ct6Qsub 


-49.6(3.4) 


21(14) 


-6(14) 


2.54 


-53.6(2.2) 


25(10) 


-2.8(1.1) 


0.41 




— asQsub 


0.6(1.1) 


-5.9(4.1) 


5.0(3.8) 


2.36 


0.04(63) 


-5.3(2.6) 


5.7(2.5) 


1.54 




Q? 


0.0555(38) 


0.057(16) 


-0.027(15) 0.18 


0.0557(17) 


0.0573(77) 


-0.0299(79) 1.50 



24' i x 32 



TABLE VI: Same as Table|3for the fit function aim 2 M + a 2 {m 2 M ) 2 + a 4 (m|f) 3 . 
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ai[GeV" 



16' x 32 
a 2 [GeV" 4 ] a 3 [GeV- 4 ] X*/doi 



Q[ 0) -0.11(20) 0.07(14) 0.02(39) 0.65 
-QiQsub -0.15(13) 0.128(90) 0.01(27) 0.20 
Qi 0) -aiQ sub -0.29(27) 0.21(19) -0.03(54) 0.13 


0.04(12) -0.110(84) 0.13(25) 2.01 
-0.369(74) 0.268(49) -0.46(16) 0.06 
-0.32(16) 0.15(11) -0.31(34) 1.05 


Qf ] -0.32(18) 0.09(13) 0.20(35) 0.09 
-Q2Qsub 0.83(13) -0.376(93) 0.01(28) 2.43 
Q^-aiQsub 0.50(21) -0.29(15) 0.21(41) 2.08 


-0.50(12) 0.259(77) 0.01(25) 0.35 
0.987(88) -0.515(56) 0.23(20) 0.52 
0.53(12) -0.276(85) 0.33(25) 1.11 


Qf ] -1.01(76) 0.32(56) 0.5(1.5) 0.38 
-a 3 Q S ub 1.22(49) -0.40(34) 0.1(1.0) 0.89 
Q!j 0) - asQsub 0.12(98) -0.01(71) 0.4(2.0) 0.56 


-0.85(45) 0.10(31) 0.59(97) 0.98 
0.86(29) -0.21(19) -0.93(66) 0.49 
0.05(58) -0.14(39) -0.2(1.2) 0.87 


Qf ] 18.7(2.4) -6.8(1.7) 4.8(5.2) 2.44 
-QsQsub -17.9(2.1) 6.1(1.4) -2.4(4.6) 2.49 
Q^-asQsub 0.5(1.1) -0.54(77) 1.9(2.1) 0.21 


19.5(1.7) -8.6(1.0) 3.1(3.8) 0.48 
-19.0(1.4) 7.99(88) -1.4(3.4) 0.38 
0.32(61) -0.46(42) 1.2(1.3) 0.61 


Qf ] 52.9(7.1) -19.5(4.9) 11(15) 2.87 
-asQsub -50.1(6.1) 16.5(4.1) -4(13) 2.58 
Qi 0) -« 6 Qsub 1.8(1.9) -2.2(1.5) 5.0(3.7) 2.35 


56.1(4.8) -24.2(3.0) 10(11) 0.32 
-54.1(4.2) 22.6(2.5) -2.5(9.7) 0.42 
1.4(1.1) -1.16(80) 5.7(2.3) 1.21 


Q <2> 0.0498(68) 0.0364(47) -0.026(15) 0.15 


0.0494(32) 0.0351(20) -0.0285(73) 0.99 



24 J x 32 

oi[GeV" 2 ] a 2 [GeV- 4 ] a 3 [GeV- 4 ] xVdof 



TABLE VII: Same as Table IV! for the fit function a\m 2 M + a2(m\ 4 ) 2 + a'j,{m 2 M ) 2 lnm|^ including a chiral logarithm term. 



m f a 0.02 


0.03 


0.04 


0.05 


0.06 


16 3 x 32 -145.0(3.1) 
24 3 x 32 -154.7(1.6) 


-135.1(2.5) 
-142.1(2.1) 


-132.6(2.3) 
-131.9(1.9) 


-120.9(1.9) 
-127.4(1.7) 


-120.2(1.8) 
-119.3(1.6) 



TABLE VIII: 



2/1 



(tt + \sd\ K + ) as a function of m/a. 
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m/a 


0.02 


0.03 


0.04 


0.05 


0.06 


bare 


<Qi)o 


-0.0329(69) 


-0.0266(48) 


-0.0242(37) 


-0.0222(29) 


-0.0176(25) 


[GeV 3 ] 




0.0508(57) 


0.0412(35) 


0.0442(27) 


0.0347(20) 


0.0345(17) 




(<9s) 


-0.006(26) 


-0.008(18) 


0.005(13) 


-0.009(10) 


0.0039(91) 




{Q4)q 


0.078(24) 


0.059(16) 


0.074(12) 


0.0475(93) 


0.0560(82) 




<Qs) 


-0.030(29) 


-0.059(18) 


-0.051(13) 


-0.054(10) 


-0.0439(88) 




(Q&) 


-0.031(62) 


-0.157(31) 


-0.109(23) 


-0.161(15) 


-0.138(12) 




(Q?)o 


1.635(30) 


2.043(33) 


2.574(42) 


2.835(43) 


3.328(49) 




(Qs) 


5.012(91) 


6.25(10) 


7.90(13) 


8.66(13) 


10.18(15) 




(Qv)o 


-0.0464(58) 


-0.0357(35) 


-0.0389(28) 


-0.0285(20) 


-0.0284(18) 




(Qio)o 


0.0372(69) 


0.0321(48) 


0.0294(38) 


0.0284(30) 


0.0237(25) 




(Qi) 2 


0.01314(15) 


0.01402(12) 


0.01487(11) 


0.015399(98) 


0.015957(90) 




<Q 2 ) 2 


0.01314(15) 


0.01402(12) 


0.01487(11) 


0.015399(98) 


0.015957(90) 




(Qr) 2 


0.4110(42) 


0.4292(34) 


0.4656(28) 


0.4863(27) 


0.5264(24) 




{Qs} 2 


1.238(13) 


1.261(11) 


1.3357(87) 


1.3639(77) 


1.4451(70) 




(Qo) a 


0.01971(23) 


0.02103(18) 


0.02231(16) 


0.02310(15) 


0.02393(13) 




(Qio) 2 


0.01971(23) 


0.02103(18) 


0.02231(16) 


0.02310(15) 


0.02393(13) 




m/a 


0.02 


0.03 


0.04 


0.05 


0.06 


renormalized 


(Qi)o 


-0.0291(68) 


-0.0234(47) 


-0.0206(37) 


-0.0191(29) 


-0.0144(25) 


at 1.3 GeV 


{Q2) 


0.0510(69) 


0.0360(43) 


0.0415(33) 


0.0291(24) 


0.0301(20) 


[GeV 3 ] 


(Q:i) 


0.004(28) 


-0.012(20) 


0.007(14) 


-0.015(11) 


0.0002(99) 






0.082(27) 


0.049(18) 


0.069(13) 


0.035(10) 


0.0460(92) 




<Qs} 


-0.026(26) 


-0.032(16) 


-0.033(12) 


-0.0253(94) 


-0.0187(82) 




{Qe) 


-0.012(48) 


-0.111(24) 


-0.071(18) 


-0.115(12) 


-0.0960(90) 




{Qr) 


0.797(17) 


1.021(18) 


1.269(21) 


1.417(21) 


1.640(23) 




{Qs)q 


3.428(69) 


4.374(73) 


5.469(86) 


6.046(87) 


7.024(94) 




<Q9> 


-0.0453(70) 


-0.0287(43) 


-0.0341(34) 


-0.0205(24) 


-0.0212(21) 




(Qio) 


0.0347(68) 


0.0306(48) 


0.0278(37) 


0.0275(29) 


0.0231(25) 




(Qi) 2 


0.01345(16) 


0.01436(13) 


0.01524(11) 


0.01578(10) 


0.016361(91) 




{Qi) 2 


0.01328(16) 


0.01417(12) 


0.01504(11) 


0.015571(99) 


0.016137(91) 




(Qs) 2 


-0.00002740(31) 


-0.00003058(27) 


-0.00003395(25) 


-0.00003677(24) 


-0.00004007(23) 




(Qa) 2 


-0.0002198(36) 


-0.0002349(30) 


-0.0002521(25) 


-0.0002652(21) 


-0.0002830(20) 




<Qs) 2 


0.0002056(37) 


0.0002196(31) 


0.0002357(25) 


0.0002483(22) 


0.0002656(20) 




<Q 6 ) 2 


0.000758(14) 


0.000789(11) 


0.0008274(91) 


0.0008517(77) 


0.0008913(70) 




{Qr) 2 


0.2045(36) 


0.2243(30) 


0.2466(25) 


0.2655(22) 


0.2897(21) 




{Qs} 2 


0.846(16) 


0.880(13) 


0.922(10) 


0.9488(86) 


0.9922(79) 




<Q 9 } 2 


0.02026(24) 


0.02161(19) 


0.02295(16) 


0.02376(15) 


0.02464(14) 




(Qio) 2 


0.02006(24) 


0.02141(19) 


0.02272(16) 


0.02353(15) 


0.02439(14) 



TABLE IX: Hadronic matrix elements {Qi) and (Qi} 2 {i = 1, • • • 10) in units of GeV 3 at each m/a on a 16 3 x 32 lattice. The 
upper half of the table lists the bare values. The lower half are those renormalized in the MS scheme at /i — 1/a and run to 
At = 1.3 GeV for N f = 3 using A ( IL = 372 MeV, which corresponds to A^L = 325 MeV. 
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0.02 


0.03 


0.04 


0.05 


0.06 


bare 


(Qi) 


-0.0235(34) 


-0.0205(37) 


-0.0217(25) 


-0.0161(22) 


-0.0196(19) 


[GeV 3 ] 


<Q2) 


0.0460(28) 


0.0457(27) 


0.0378(19) 


0.0351(16) 


0.0345(11) 




(<93> 


0.013(12) 


0.021(13) 


0.0010(98) 


0.0116(81) 


0.0003(68) 




{Q4)q 


0.082(11) 


0.087(12) 


0.0600(91) 


0.0627(74) 


0.0544(58) 




(Qs) n 


-0.027(14) 


-0.044(13) 


-0.0515(97) 


-0.0393(82) 


-0.0442(65) 




(<3e) 


-0.105(26) 


-0.183(24) 


-0.167(17) 


-0.136(14) 


-0.117(11) 




(Qr) 


1.697(15) 


2.157(29) 


2.563(35) 


2.990(40) 


3.295(44) 




(Qs) n 


5.211(44) 


6.584(85) 


7.84(11) 


9.13(12) 


10.08(13) 






-0.0417(28) 


-0.0412(28) 


-0.0324(20) 


-0.0299(17) 


-0.0295(12) 




(Qio) 


0.0278(34) 


0.0250(38) 


0.0267(26) 


0.0212(22) 


0.0246(19) 




(<?1>2 


0.013154(43) 


0.014163(52) 


0.014781(48) 


0.015335(45) 


0.015853(43) 




<<52> 2 


0.013154(43) 


0.014163(52) 


0.014781(48) 


0.015335(45) 


0.015853(43) 






0.3996(15) 


0.4222(18) 


0.4559(15) 


0.4900(14) 


0.5184(13) 




<Qs) 2 


1.2119(48) 


1.2444(55) 


1.3128(45) 


1.3783(41) 


1.4271(40) 




<Qo) 2 


0.019730(65) 


0.021244(78) 


0.022172(72) 


0.023003(67) 


0.023779(65) 




(Qio) 2 


0.019730(65) 


0.021244(78) 


0.022172(72) 


0.023003(67) 


0.023779(65) 




m/a 


0.02 


0.03 


0.04 


0.05 


0.06 


renormalized 


(Qi)o 


-0.0203(34) 


-0.0173(36) 


-0.0182(25) 


-0.0130(21) 


-0.0163(19) 


at 1.3 GeV 




0.0433(33) 


0.0401(34) 


0.0322(23) 


0.0309(19) 


0.0310(14) 


[GeV 3 ] 


<Qs) 


0.015(14) 


0.017(14) 


-0.004(11) 


0.0084(88) 


-0.0014(73) 




(04) o 


0.078(13) 


0.076(14) 


0.048(10) 


0.0535(84) 


0.0466(65) 




<Qs) 


-0.008(12) 


-0.011(11) 


-0.0214(90) 


-0.0144(74) 


-0.0231(57) 






-0.072(20) 


-0.132(19) 


-0.120(14) 


-0.095(11) 


-0.0790(84) 




(Qr)o 


0.8415(80) 


1.072(15) 


1.271(17) 


1.488(19) 


1.637(21) 




(<3s) 


3.631(33) 


4.566(60) 


5.434(70) 


6.352(79) 


7.002(86) 






-0.0376(33) 


-0.0338(34) 


-0.0247(24) 


-0.0233(20) 


-0.0231(14) 




(Qio)o 


0.0259(34) 


0.0234(37) 


0.0256(26) 


0.0205(21) 


0.0239(19) 




<<3l>2 


0.013469(44) 


0.014499(53) 


0.015140(49) 


0.015717(46) 


0.016253(44) 




(<52> 2 


0.013295(44) 


0.014317(53) 


0.014944(49) 


0.015507(45) 


0.016031(43) 




<Q3) 2 


-0.00002694(11) 


-0.00003018(14) 


-0.00003331(13) 


-0.00003662(13) 


-0.00003960(13) 




<Q 4 ) 2 


-0.0002180(15) 


-0.0002301(17) 


-0.0002476(15) 


-0.0002660(13) 


-0.0002807(13) 




{Qs) 2 


0.0002037(15) 


0.0002142(17) 


0.0002312(15) 


0.0002492(13) 


0.0002634(13) 




(Q 6 ) 2 


0.0007562(57) 


0.0007728(63) 


0.0008144(54) 


0.0008575(48) 


0.0008866(46) 




<Q 7 ) 2 


0.2010(15) 


0.2180(17) 


0.2409(15) 


0.2658(14) 


0.2866(14) 




<Qs) 2 


0.8440(64) 


0.8618(70) 


0.9078(60) 


0.9552(53) 


0.9872(51) 




<Qs)) 2 


0.020281(66) 


0.021830(80) 


0.022796(74) 


0.023666(69) 


0.024474(66) 




(Qio) 2 


0.020083(66) 


0.021623(80) 


0.022575(73) 


0.023431(68) 


0.024228(66) 



TABLE X: Same as Table llXl for the 24 3 x 32 lattice. 
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g* = l/o 


q* — ir/a 


(Qi) 


-0.0203(34) 


-0.0152(33) 


{Q2) 


0.0433(33) 


0.0424(33) 


(Qs) 


0.015(14) 


0.019(14) 


<Q4} 


0.078(13) 


0.076(13) 


<Qs) 


-0.008(12) 


-0.005(12) 


<Q6> 


-0.072(20) 


-0.050(16) 


<Q 7 ) 


0.8415(80) 


0.7986(77) 


{Qs) 


3.631(33) 


2.873(26) 


(Q9) 


-0.0376(33) 


-0.0317(33) 


(Qio)o 


0.0259(34) 


0.0257(34) 


<Ql>2 


0.013469(44) 


0.014314(46) 


<Q 2 } 2 


0.013295(44) 


0.013760(45) 


<Qt} 2 


0.2010(15) 


0.1912(14) 


<Q 8 ) 2 


0.8440(64) 


0.6678(51) 


<Q 9 } 2 


0.020281(66) 


0.021754(70) 


<QlO> 2 


0.020083(66) 


0.021149(69) 



TABLE XI: Renormalized hadronic matrix elements at fi = 1.3 GeV in units of GeV 3 from different matching point q* = 1/a 
(left column) and n/a (right column). Values are taken at m/a = 0.02 on a 24 x 32 lattice. 





quadratic 


X 2 /dof 


chiral log. 


X 2 /dof 


(Qi) 


-0.034(20) 


0.14 


-0.035(36) 


0.14 


<Q 2 ) 


0.070(19) 


3.03 


0.083(34) 


3.05 


(Q 3 > 


0.033(82) 


0.91 


0.06(15) 


0.94 


(Qi) 


0.131(78) 


1.68 


0.17(14) 


1.71 


<Q 5 ) 


-0.008(72) 


0.03 


0.03(13) 


0.02 


<Q 6 ) 


0.08(12) 


2.64 


0.20(21) 


2.63 


(Qs) (4pts.) 


-0.04(17) 


4.32 


-0.02(31) 


4.38 


(Qr) 


0.247(78) 


1.78 


0.11(15) 


1.69 


(Qs) 


1.07(32) 


2.87 


0.48(60) 


2.77 


(Q9} 


-0.067(19) 


3.32 


-0.082(35) 


3.35 


(Qio)o 


0.037(21) 


0.17 


0.037(37) 


0.17 


<Ql>2 


0.01102(54) 


0.34 


0.00990(98) 


0.29 


<Q 2 ) 2 


0.01087(54) 


0.33 


0.00975(97) 


0.28 


<Q3> 2 


-0.0000203(12) 


0.49 


-0.0000195(21) 


0.46 


<Q 4 ) 2 


-0.000188(12) 


0.33 


-0.000187(22) 


0.37 


<Q 5 ) 2 


0.000177(12) 


0.33 


0.000179(23) 


0.34 


<Qe) 2 


0.000694(46) 


0.31 


0.000694(83) 


0.31 


{Qr} 2 


0.164(12) 


0.36 


0.167(23) 


0.38 


(Qs} 2 


0.776(51) 


0.30 


0.775(92) 


0.31 




0.01660(81) 


0.34 


0.0149(15) 


0.29 


(Qw} 2 


0.01642(81) 


0.33 


0.0147(15) 


0.29 



TABLE XII: Hadronic matirx elements in units of GeV 3 in the chiral limit m 2 M — > on a 16 3 x 32 lattice. The column 
named "quadratic" , "chilal log." or "(quadratic) + breaking" corresponds to three types of fit forms described in the text. 
Chiral extrapolations are made using data at all m/a = 0.02-0.06 (5 points) except for an alternative extraporation of {Qe) 
exclluding the point at m/a = 0.02(4 points). 
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quadratic 


X 2 /dof 


chiral log. 


X 2 /dof 


(Qi)o 


-0.031(12) 


0.98 


-0.039(21) 


1.04 


<Q 2 ) 


0.066(11) 


0.36 


0.082(20) 


0.43 


<Q3) 


0.032(48) 


0.83 


0.044(86) 


0.82 


{Qi) 


0.124(45) 


0.69 


0.155(82) 


0.69 


(Q5> 


0.001(40) 


0.44 


0.003(73) 


0.43 


(Qe) 


0.014(66) 


1.53 


0.16(12) 


1.20 


(Q 6 ) (4pts.; 


-0.19(13) 


0.10 


-0.18(25) 


0.10 


{Qr) 


0.252(53) 


0.43 


0.07(11) 


0.45 


(Qs> 


1.23(22) 


0.27 


0.58(43) 


0.35 


{Qs)o 


-0.063(11) 


0.45 


-0.082(20) 


0.53 


(Qio} 


0.034(12) 


1.09 


0.039(22) 


1.14 


<Ql> 2 


0.01104(19) 


4.55 


0.00979(36) 


2.99 


<Q2>2 


0.01089(19) 


4.79 


0.00964(35) 


3.16 


(Qb) 2 


-0.00001942(50) 


0.25 


-0.00001832(96) 


0.18 


<Q4> 2 


-0.0001848(59) 


1.28 


-0.000187(11) 


1.25 


<Q 5 ) 2 


0.0001737(61) 


1.67 


0.000179(11) 


1.57 


<Qe} 2 


0.000691(22) 


2.08 


0.000708(42) 


1.98 


<Qr) 2 


0.1580(60) 


1.27 


0.163(11) 


1.18 


<Qs} 2 


0.772(25) 


2.09 


0.792(47) 


1.99 


(Q9) 2 


0.01663(28) 


4.48 


0.01476(54) 


2.94 


(Qio} 2 


0.01646(28) 


4.66 


0.01458(53) 


3.07 



TABLE XIII: Same as Table |XIT| for the 24 3 x 32 lattice. 





16 3 x 32 




24 3 x 32 






quadratic 


chiral log. 


quadratic 


chiral log. 


S (V2) 


8.3(5.0) 


8.6(8.9) 


7.7(2.9) 


9.6(5.2) 


R (l/2) 


3.43(95) 


4.1(1.7) 


3.23(55) 


4.04(98) 


R (l/2) 


2.7(6.7) 


5(12) 


2.6(3.9) 


3.6(7.1) 


R (l/2) 


3.6(2.1) 


4.5(3.8) 


3.4(1.2) 


4.3(2.3) 


R (l/2) 
^5 


0.04(40) 


-0.15(71) 


0.01(22) 


-0.02(41) 


B (l/2) 


-0.14(22) 


-0.38(38) 


-0.03(12) 


-0.29(22) 


^ 1/2) (4 P ts.; 


0.07(31) 


0.03(58) 


0.35(25) 


0.34(47) 




0.49(15) 


0.22(29) 


0.50(10) 


0.14(21) 


R (l/2) 


0.73(22) 


0.32(41) 


0.83(15) 


0.39(29) 


R (l/2) 


5.5(1.6) 


6.8(2.8) 


5.19(92) 


6.7(1.7) 


r(1/2) 


3.0(1.7) 


3.0(3.0) 


2.78(98) 


3.2(1.8) 


fl (3/2) 


0.480(24) 


0.431(43) 


0.4809(82) 


0.426(16) 


fl (3/2) 


0.473(23) 


0.425(42) 


0.4745(81) 


0.420(15) 


S (3/2) 


0.640(49) 


0.651(88) 


0.616(23) 


0.634(44) 


r(3/2) 


0.924(61) 


0.92(11) 


0.920(30) 


0.944(55) 


r(3/2) 


0.482(24) 


0.433(43) 


0.4830(82) 


0.429(16) 


r(3/2) 
D 10 


0.477(24) 


0.428(43) 


0.4779(82) 


0.423(15) 



TABLE XIV: B parameters in the chiral limit with the chiral logarithm fit. 
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16 3 x 32 


24 3 x 32 






ReA [10~ 8 GeV] ReA 2 [10" 


*GeV] oj- 1 


ReA [10- 8 GeV] ReA 2 [10" 8 GeV] 




A (4) 

MS 


= 325 MeV 
















0.02 


13.1(1.4) 


1.867(22) 


7.01(78) 


10.80(69) 


1.8689(62) 


5.78(37) 




0.03 


9.45(84) 


1.992(17) 


4.75(43) 


9.90(69) 


2.0129(74) 


4.92(35) 




0.04 


10.42(68) 


2.114(15) 


4.93(33) 


8.26(45) 


2.1006(68) 


3.93(22) 




0.05 


7.66(49) 


2.188(14) 


3.50(22) 


7.61(38) 


2.1792(64) 


3.49(17) 




0.06 


7.52(40) 


2.267(13) 


3.32(17) 


7.86(28) 


2.2527(61) 


3.49(12) 


A (4) 

MS 


= 215 MeV 
















0.02 


12.5(1.4) 


1.911(23) 


6.55(72) 


10.40(66) 


1.9130(63) 


5.43(34) 




0.03 


9.12(81) 


2.039(18) 


4.47(40) 


9.59(66) 


2.0602(76) 


4.66(32) 




0.04 


10.04(64) 


2.164(16) 


4.64(31) 


8.01(43) 


2.1500(70) 


3.72(20) 




0.05 


7.41(47) 


2.240(14) 


3.31(21) 


7.38(36) 


2.2306(65) 


3.31(16) 




0.06 


7.30(38) 


2.321(13) 


3.15(16) 


7.60(26) 


2.3058(63) 


3.29(12) 


MR 


= 435 MeV 
















0.02 


13.7(1.5) 


1.821(22) 


7.52(84) 


11.20(72) 


1.8228(60) 


6.14(40) 




0.03 


9.78(89) 


1.943(17) 


5.03(46) 


10.18(73) 


1.9635(72) 


5.19(38) 




0.04 


10.80(71) 


2.062(15) 


5.24(35) 


8.50(48) 


2.0489(67) 


4.15(23) 




0.05 


7.87(51) 


2.134(14) 


3.69(24) 


7.82(40) 


2.1254(62) 


3.68(19) 




0.06 


7.74(42) 


2.211(12) 


3.50(19) 


8.11(29) 


2.1970(60) 


3.69(13) 



TABLE XV: Values of ReA , RcA 2 and lo 1 obtained at each rata for both lattice sizes, with A< 4 1=325 MeV, 215 MeV and 

' MS 

435 MeV. 



24 3 x 32 




£o 


£i 


6 


6 


X 2 /dof 


quadratic Re^lo [10~ 


j GeV] 


16.5(2.2) 


-27.7(9.2) 




21.8(8.8) 


0.34 


Re^ 2 [10~ 


^GeV] 


1.531(26) 


1.62(12) 




-0.86(13) 


4.91 


CO' 1 




9.5(1.1) 


-18.2(4.6) 




13.7(4.3) 


0.13 


chiral log. RcA) [10" 


3 GeV] 


20.7(4.0) 


-11.4(2.8) 


20.1(8.3) 




0.50 


Re^ 2 [10" 


3 GeV] 


1.353(50) 


0.977(31) 


-0.82(11) 




3.25 


LO- 1 




12.3(2.0) 


-8.0(1.5) 


12.9(4.1) 




0.26 



TABLE XVI: Fit parameters for ReAo, Re^b and to 1 with £ + $,im 2 M + S,z{rn 2 M ) 2 (quadratic fit) and £o +£i m M +£,2rn 2 M \nm 2 M 
(chiral logarithm fit). Results on a 24 3 x 32 lattice with A^L = 325 MeV are shown. 
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16 3 x 32 24 3 x 32 


p(l/2) p(3/2) £ , /e [10 - 4] 


p(l/2) p(3/2) £ , /e [10 -4j 


A^L = 325 MeV 

MR 


0.02 0.1(1.2) 4.93(11) -6.3(1.5) 
0.03 2.97(61) 5.084(88) -2.74(78) 
0.04 2.19(44) 5.291(70) -4.03(56) 
0.05 3.65(29) 5.416(59) -2.30(37) 
0.06 3.42(22) 5.657(53) -2.90(28) 


1.69(50) 4.923(45) -4.21(64) 
3.36(47) 4.944(49) -2.06(60) 
3.50(33) 5.200(41) -2.21(41) 
3.19(26) 5.470(37) -2.96(33) 
3.01(19) 5.632(35) -3.41(24) 


AH. = 215 MeV 

MR 


0.02 0.06(94) 3.713(87) -4.7(1.2) 
0.03 2.38(50) 3.815(70) -1.86(63) 
0.04 1.74(36) 3.962(56) -2.89(45) 
0.05 2.93(24) 4.049(47) -1.45(30) 
0.06 2.75(18) 4.228(42) -1.92(22) 


1.34(41) 3.707(36) -3.07(52) 
2.70(38) 3.701(39) -1.31(49) 
2.81(26) 3.892(33) -1.40(33) 
2.56(21) 4.094(27) -1.99(27) 
2.41(16) 4.211(28) -2.34(19) 


A^L = 435 MeV 

MR 


0.02 0.1(1.4) 6.16(13) -7.8(1.8) 
0.03 3.63(75) 6.36(11) -3.56(95) 
0.04 2.67(54) 6.629(84) -5.15(68) 
0.05 4.43(36) 6.790(71) -3.06(46) 
0.06 4.16(27) 7.091(64) -3.82(34) 


2.05(61) 6.150(54) -5.33(78) 
4.09(58) 6.197(58) -2.74(73) 
4.26(40) 6.518(50) -2.93(50) 
3.88(32) 6.853(44) -3.87(40) 
3.65(24) 7.059(42) -4.43(29) 



TABLE XVII: Values of pt 1 ^ pO/2) and £ y £ at each mja for both lattice volumeS) with A^.^325 MeV, 215 MeV and 435 
MeV. 
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FIG. 1: (Left) Anomalous quark mass ms, as a function of N5 in the m/a — > limit for the RG-improved action. Filled (open) 
circles represent data at ((3, M) = (2.2,1.7) on a 16 3 (12 3 ) x 24 lattice. Filled squares are those at (/3, M) = (2.6,1.8) on a 
16 3 x 32 lattice. For the latter, data at four larger N5 are used for fits with the functions ae~^ N5 (dotted line) and c + ae^ 5 
(solid line). (Right) Same for the plaquette action at (/?, M) =(5.65, 1.7) and (6.0, 1.8). 
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t 



(a) figure-eight 

K 



(b) eye 




K 



<ji + IQ ( ?IKS 




<7T + IQ ( ?IK + > 



(c) annihilation 



(type 1) K 



(type 2) K u 




-^-<0IQ ( °k°> 
dm s 



(d) subtraction 



K 





<7T + IQ sub IK + > 



(e) two-point 
K 




<0IQ sub IK°> 



FIG. 2: Types of contractions needed for our calculation. Solid lines represent quark propagators on a background gauge field. 
Crosses represent points where meson sources are placed, while filled squares denote four quark operators or the subtraction 



operator, (a) "figure-eight", (b) "eye" which contributes only for matrix elements of Q 



(0) 



" annihilation" with a quark mass 



derivative in the external line (type f) or in the quark loop (type 2), (d) "subtraction" and (e) "two-point". 
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FIG. 3: Time dependence of the propagator ratio defined by 111! .2211 for (upper) and Q 6 (lower) for mja = 0.03. Left 
and right columns are for the lattice size 16 3 x 32 and 24 3 x 32 respectively. 




FIG. 4: Time dependence of the propagator ratio denned by lLi.1.251 to calculate the parameter a 2 x 02 (upper) and a 2 x 
(lower) at m/a = 0.03. Left and right columns are for the lattice size 16 3 x 32 and 24 3 x 32 respectively. 
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FIG. 5: Effect of subtractions illustrated for (upper) and Q 1 ^ 1 (lower) as a function of m\ { . The original matrix element 



(0) 



)(<>) 



K + J (circles) and the subtraction term — cti (n + jQ su b| ^ + ) (diamonds) are added to obtain the physical matrix 

element (squares). Values are multiplied with a factor \/2f n (m 2 K — m'i)/ \ A* \ 0) (0 \ A± \ so that the vertical axis has 
dimension [GeV 3 ] . Left and right columns are for the lattice sizes 16 3 x 32 and 24 3 x 32 respectively. 
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FIG. 



Ratio of matrix elements 



X 



(I) 



/(tt+ \A4\O) (0\A 4 \K + ) x m 2 M a 2 as a function of m|f[GeV 2 ] for 



1, 2, 3, 5, 6 (J = 0) and i = 1 (7 = 2) from top to bottom. Left and right columns are for the lattice sizes 16 x 32 and 24 x 32 



respectively. Solid lines represent the chiral extrapolation to m 2 M 
with a cubic function as described in the text. 



with a quadratic function of m M , while dashed lines are 



0.00 



-0.02 



-0.04 



-0.06 



-0.08 



<Q,>o 



□ 16 x32 
■ 24 3 x32 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
m M 2 [GeV 2 ] 



-, 0.12 




<Q 2 > 


10 






0.08 




□ 16 3 x32 






»24 3 x32 h 


- 0.06 






0.04 












0.02 






- 0.00 







0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

m M 2 [GeV 2 ] 



0.2 - 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
m M 2 [GeV 2 ] 



0.15 
0.10 
0.05 
0.00 
-0.05 
-0.10 
-0.15 



<Q 5 >o 



□ 16 x32 
■ 24 3 x32 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
m M 2 [GeV 2 ] 



0.2 
0.1 
0.0 
-0.1 
-0.2 
-0.3 
-0.4 
-0.5 





<Q6>0 












□ 1 6 3 x32 . 






■ 24 3 x32 









0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
m M 2 [GeV 2 ] 



44 



0.08 




<Qio>o 




0.00 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 



< [GeV 



m M [Gev] 



FIG. 7: Physical hadronic matrix elements (Q;} for i = 1,2,3,4,5,6,9,10 as a function of m 2 M from top to bottom. These 
matrix elements involve subtractions of unphysical effects. Open and filled symbols are from the spatial volume V = 16 3 and 
24 3 respectively. Chiral extrapolations with a quadratic polynomial is shown by solid (V = 24 3 ) and dashed (V = 16 3 ) lines. 
Fit error in the chiral limit is added for the former. 
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FIG. 9: Physical hadronic matrix elements, (Qi} 2 and (Q7,s} 2 as a function of m 2 M . The organization of each panel is the same 
as that in Fig. Q 
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FIG. 10: ReA2 (left) and ReAo (right) in units of GeV as a function of m 2 M . For chiral extrapolation, quadratic (solid) and 
chiral logarithm (dashed) forms are used. For the formar, fit errors are shown in the chiral limit. Filled and open symbols are 
for the spatial volume 24 3 and 16 3 respectively. 
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FIG. 11: Breakdown of ReA2 (left) and KeAo (right) into contributuions from the opertors Qi(i = 1, • • • , 10) at m;o = 0.03. 
Data points placed on horizontal lines show total values and errors. The solid and dashed lines are for the spatial volume 24 3 
and 16 3 respectively. 
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FIG. 12: KeAo/ReA2 as a function of m 2 M . For chiral extrapolation, quadratic form (solid line with the fit error at m 2 M = 0) 
and chiral logarithm form (dashed line) are used. Open and filled symbols are for the spatial volume 16 3 and 24 3 respectively. 
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FIG. 13: P (3/2) (left) and P (1/2) (right) as a ftmction of m M . Open and filled symbols are for the spatial volume 16 and 24 
respectively. 
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FIG. 14: Breakdown of P (3/2) (left) and P (1/2) (right) into contributuions from the opertors Q;(i = 3, • • • , 10) at m f a = 0.03. 
Data points placed on horizontal lines show total values and errors. The solid and dashed lines are for the spatial volume 24 3 
and 16 a respectively. 
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FIG. 15: e'/e as a function of m? M . Open and filled symbols are for the spatial volume 16 3 and 24 3 respectively. Experimental 
values quoted in (ll.4t are also shown. 



